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1 1  ABSTRACT  no  J  f 

In  multichannel  identification  and  detection  (or  model-based  multichannel  detection) 
problem*  the  parameters  of  a  model  are  identified  from  the  observed  channel  process, 
and  the  identified  model  is  used  to  facilitate  the  detection  of  a  desired  signal  in 
the  observed  process*  A  model-based  multichannel  detection  algorithm  was  developed  in 
the  context  of  an  innovations-based  detection  algorithm  (IBDA)  formulation  for  sur¬ 
veillance  radar  system  applications.  The  state  space  model  class  was  adopted  to  model 
the  vector  channel  process  because  it  is  more  general  than  the  time  series  model  class 
used  in  most  analyses  to  date.  An  IBDA  methodology  was  developed  based  on  an  algorithm 
which  uses  output  data  directly  and  offers  computational  and  performance  advantages 
over  alternative  techniques*  A  computer  simulation  was  developed  to  validate  the 
methodology  and  the  algorithm,  and  to  carry  out  performance  assessments.  Simulation 
results  indicate  that  the  algorithm  is  capable  of  discriminating  between  the  null 
hypothesis  (clutter  plus  noise)  and  the  alternative  hypothesis  (signal  plus  clutter 
plu*  noise).  In  summary,  the  applicability  of  the  approach  to  radar  system  problems  has 
been  established. 
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1 . 0  INTRODUCTION 


In  multichannel  identification  problems  the  outputs  of 
multiple  channels  (or  sensors)  are  available,  and  it  is  desired  to 
identify  the  parameters  of  an  analytical  model  to  represent  the 
phenomena  being  observed  via  the  channel  outputs.  Similarly,  in 
multichannel  detection  problems  the  outputs  of  multiple  channels 
are  available,  and  it  is  desired  to  determine  the  presence  (or 
absence)  of  a  desired  signal  component  in  the  channel  data.  In 
the  combined  problem  of  multichannel  identification  and  detection 
a  model  is  estimated  for  the  phenomena  being  observed  via  the 
channel  outputs,  and  the  identified  model  is  used  to  facilitate 
the  detection  of  a  desired  signal  in  the  channel  output  data. 
Multichannel  identification  and  detection  is  thus  referred  to  also 
as  model-based  multichannel  detection.  In  all  of  these  problems 
the  channel  data  is  available  simultaneously  over  many  channels  of 
the  same  type,  or  over  many  distinct  channels  (each  channel 
corresponding  to  a  different  sensor  type) . 

This  report  is  a  summary  of  the  work  carried  out  in  Phase  I 
of  this  program.  Specifically,  the  development  of  state  space 
algorithms  for  model-based  multichannel  detection  in  the  context 
of  surveillance  radar  system  applications  is  addressed.  In 
surveillance  radar  systems  (radar  arrays)  the  channels  correspond 
to  separate  antenna  apertures  (or  elements  of  a  single  aperture 
array)  .  The  desired  signal  may  or  may  not  be  present  in  the 
channel  output  data  at  any  given  time.  The  data  in  each  channel 
generally  includes  noise  (broadband  interference)  as  well  as 
"clutter"  (narrowband  interference),  with  low  s ignal-to-clutter 
ratio  and,  possibly,  low  signal-to-noise  ratio  also.  Model-based 
detection  methods  must  discriminate  between  the  condition  of 
target  embedded  in  clutter  and  noise,  and  the  condition  of  clutter 
and  noise  only. 
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Figure  1-1  illustrates  a  radar  array  system  consisting  of 
multiple  subarrays  or  array  elements.  The  output  of  each  subarray 
(or  each  individual  array  element)  is  a  complex-valued,  scalar, 
digital  sequence,  denoted  as  (Xj(n)}.  The  collection  of  the  J  scalar 

sequences  is  arranged  into  a  »)-dimensional  vector,  {&(n)},  which  is 
input  to  a  multichannel  processor  (not  shown  in  the  figure) . 


Channel  No.  1 


Channel  No.  J 


Figure  1-1.  -Radar  array  with  J  subarrays  or  individual  elements. 

In  Phase  I  the  multivariate  (multiple  input,  multiple  output) 
state  space  model  class  was  adopted  to  represent  the  multichannel 
radar  data,  and  new  system  identification  techniques  were  applied 
to  estimate  the  model  parameters.  The  modeling  of  the  complex¬ 
valued  pre-processed  radar  signals  for  multichannel  detection 
using  the  state  space  model  class  is  one  of  the  contributions  of 
this  work.  State  space  models  have  been  used  in  the  context  of 
target  tracking  (where  the  detected  radar  signal  is’  processed 
further  to  estimate  a  trajectory)  and  for  the  determination  of 
weights  in  antenna  array  sidelobe  canceling  and  related  problems, 
but  not  for  multichannel  detection.  Model-based  detection  has 
been  carried  out  using  the  more-restricted  time  series  models 
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(Michels,  1991;  Metford  and  Haykin,  1985),  which  are  included 
within  the  class  of  state  space  models  and  can  be  represented  as 
such. 


The  state  space  identification  algorithm  on  which  the 
methodology  developed  in  this  program  is  based  has  several  unique 
features.  Foremost  among  these,  the  algorithm  operates  on  output 
data  directly  to  generate  estimates  of  the  parameters  of  a  state 
space  model  (without  computing  output  correlation  matrices) .  This 
feature  of  the  algorithm  results  in  reduced  dynamic  range 
requirements  in  comparison  with  state  space  algorithms  that 
operate  on  correlation  matrices.  The  algorithm  belongs  to  the 
class  referred  to  as  subspace  methods  because  the  fundamental 
operation  of  the  algorithm  is  to  decompose  the  vector  space 
spanned  by  the  channel  output  data  into  signal  and  noise 
subspaces.  Implementation  of  this  fundamental  operation  is 
carried  out  using  the  QR  decomposition  and  the  singular  value 
decomposition  (SVD) ,  which  are  stable  numerical  techniques.  This 
identification  algorithm  is  new;  it  is  scheduled  to  appear  in  the 
open  literature  late  this  year  (Van  Overschee  and  De  Moor,  1993) . 

An  important  distinction  in  the  context  of  radar  system 
applications  is  that  the  vector  random  processes  which  represent 
the  channel  data  are  complex-valued  processes  in  most  cases.  Most 
time  series  techniques  and  models  have  been  formulated  for  complex 
as  well  as  real  processes.  The  same,  however,  cannot  be  said 
about  state-space  techniques;  state-space  methods  and  results 
available  in  the  literature  have  been  defined  almost  exclusively 
for  the  case  of  real-valued  processes,  including  the  algorithm  of 
Van  Overschee  and  De  Moor  (1993) .  In  Phase  I  the  Van  Overschee-De 
Moor  algorithm  was  extended  to  the  case  of  complex-valued 
processes,  which  is  the  formulation  presented  in  this  report. 
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A  computer  simulation  was  generated  as  part  of  this  program 
to  validate  the  methodology  and  the  algorithms,  and  to  carry  out 
simulation-based  analyses.  This  software  was  exercised  with 
simulated  multichannel  data  generated  at  RL,  and  the  modeling  and 
identification  results  compare  favorably  with  the  results  obtained 
at  RL  using  auto-regressive  models. 

In  summary,  the  analytical  and  simulation  results  obtained  in 
this  program  indicate  that  the  SSC  algorithm  and  methodology  for 
model-based  multichannel  detection  has  the  potential  to  result  in 
significant  advances  for  radar  system  applications. 

i .  i  aatatloa 

Vector  variables  are  denoted  by  underscored  lower-case 
letters  (including  Greek  letters) .  Matrices  are  denoted  by  upper¬ 
case  letters  (including  Greek  letters) .  Some  scalars  (such  as  the 
order  of  the  state  variable  model)  are  denoted  also  by  upper-case 
letters.  Vector  spaces  are  denoted  by  upper-case  script  letters, 
such  as  The  expectation  operator  is  denoted  as  E[*];  superscript 

T  and  H  are  used  to  denote  the  matrix  and  vector  transpose  and  the 
Hermitian  transpose  operators,  respectively;  and  an  asterisk  (*) 
denotes  the  complex  conjugate  operator.  1^  denotes  an  M- 
dimensional  identity  matrix,  O^j  denotes  an  NxJ  null  (zero) 
matrix,  0^  denotes  an  M-dirnensional  (square)  null  matrix,  and 
denotes  an  M-dimensional  zero  vector.  |A|  denotes  the  determinant 
of  matrix  A;  A'1  denotes  the  inverse  of  matrix  A;  At  denotes  the 
pseudoinverse  of  A;  range (A)  denotes  the  range  (column  space)  of 
A;  rank  (A)  denotes  the  rank  of  A;  A(i,j)  and  Rjj  are  both  used  to 
denote  the  (i ,j)t h  element  of  matrix  A;  and  dim ('l'')  denotes  the 
dimension  of  vector  space  V.  A  caret  (A)  over  a  variable  denotes 

an  estimate  of  the  variable,  a  bar  (— )  over  a  variable  is  used  to 
represent  the  mean  of  the  variable,  and  In (a)  denotes  the  natural 
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logarithm  of  a.  The  symbol  1  denotes  "is  orthogonal  to;"  n 
denotes  intersection  of  two  vector  spaces;  ©  denotes  the  direct 
sum  of  vector  spaces;  V  denotes  "for  all;"  and  €  denotes  "is  an 
element  of." 

Where  possible,  the  symbols  used  herein  to  represent 
variables  match  the  symbols  used  by  Michels  (1991)  to  facilitate 
enhancing  the  software  available  at  Rome  Laboratory  (RL)  with  the 
techniques  developed  in  this  program.  This  philosophy  forces  the 
use  of  non-standard  symbols  to  represent  the  parameters  of  a  state 
variable  model.  Of  course,  notational  convention  should  not  be  a 
major  issue  provided  ail  symbols  are  defined  appropriately. 
However,  it  is  important  to  mention  this  point  in  order  to  avoid 
possible  confusion  on  the  part  of  the  reader. 

l .  2  fttpart  QYirYltw 

An  introduction  to  the  model-based  multichannel  detection 
problem  is  presented  in  Section  2.0.  This  section  includes  also 
the  definition  of  the  state  space  model  class  and  several  related 
concepts,  including  the  backward  model  associated  with  a  forward 
model,  and  the  innovations  representation  for  a  random  process. 
The  parameter  identification  algorithm  is  presented  in  Section 
3.0,  and  the  algorithm  proof  provided  differs  significantly  from 
the  proof  given  by  Van  Overschee  and  De  Mtor  (1993) .  In  fact,  the 
algorithm  proof  given  here  is  simpler  and  easier  to  follow,  As 
mentioned  earlier,  this  algorithm  is  the  backbone  of  the 
Scientific  Studies  Corporation  (SSC)  multichannel  detection 
approach.  Kalman  filtering  of  the  channel  data  to  generate  the 
innovations  sequence  is  discussed  in  Section  4.0.  The  innovations 
sequence  is  fed  to  a  likelihood  ratio  detector  which  generates  the 
detection  decision,  as  described  in  Section  5.0,  A  discussion  of 
the  software  generated  in  the  program  is  presented  in  Section  6.0, 
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along  with  several  simulation  results.  Section  7.0  includes  the 
main  conclusions  and  recommendations  borne  out  of  this  Phase  I. 
Appendix  A  presents  a  methodology  for  generating  the  state  space 
representation  of  three  conventional  time  series  models  (moving- 
average,  auto-regressive,  and  auto-regressive  moving-average). 
Appendix  B  presents  the  quotient  singular  value  decomposition 
(QSVD)  for  matrix  pairs,  as  required  in  Section  3.0. 


2 . 0  MODEL-BASED  MULTICHANNEL  DETECTION 


The  model-based  approach  to  multichannel  detection  involves 
processing  the  channel  data  to  identify  the  parameters  of  a  model 
for  the  multichannel  system,  and  determination  of  a  detection 
decision  utilizing  the  identified  parameters  to  filter  the  channel 
data.  Model  parameters  can  be  identified  on-line,  as  the. channel 
data  is  received  and  processed.  Alternatively,  the  model 
parameters  can  be  identified  off-line  for  various  conditions  and 
stored  in  the  processor  memory  to  be  accessed  in  real-time  as 
required. 

There  are  two  general  classes  of  linear  parametric  models  for 
vector  random  processes:  time  series  models  and  state  space 
models.  Time  series  models  include  moving-average  (MA)  models, 
auto-regressive  (AR)  models,  and  auto-regressive  moving-average 
(ARMA)  models.  State  space  models  are  more  general  than  time 
series  models;  in  fact,  MA,  AR,  and  ARMA  models  can^be  represented 
by  state  space  models  (Appendix  A)  .  In  the  state  space 
literature,  the  determination  of  the  model  parameters  based  on 
output  data  (and,  sometimes,  input  data  also)  is  referred  to  as  a 
stochastic  identification  or  a  stochastic  realization  problem. 

Time  series  models  have  been  applied  to  the  multichannel 
detection  problem,  and  the  performance  results  obtained  provide 
encouragement  for  further  research  (see,  for  example,  Michels, 
1991,  and  the  references  therein).  The  results  obtained  by 
Michels  (1991)  assume  that  the  multichannel  output  process  can  be 
modeled  as  a  vector  AR  process.  Given  the  generality  of  state- 
space  models  and  the  wealth  of  results  available  in  the  state- 
space  literature,  the  state  space  model  class  was  selected  in 
Phase  I  to  represent  the  multichannel  signals  in  the  model-based 
multichannel  detection  problem  for  radar  systems. 
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In  the  case  of  time  series  models,  two  types  of  model 
parameter  estimation  algorithms  have  been  established  in  the 
literature:  (a)  algorithms  which  operate  on  channel  output 
correlation  matrices,  such  as  the  extended  Levinson  algorithm 
(Anderson  and  Moore,  1979),  and  (b)  algorithms  which  operate  on 
the  channel  output  data  directly  (without  the  need  to  compute 
channel  output  correlation  matrices) ,  such  as  the  Levinson- 
Wiggins-Robinson  algorithm  (Wiggins  and  Robinson,  1965)  and  the 
Strand-Nuttall  algorithm  (Strand,  1977;  Nuttall,  1976). 

In  the  case  of  state-space  models,  most  of  the  existing 
algorithms  operate  on  channel  output  correlation  matrices,  such  as 
the  stochastic  realization  approach  developed  by  Akaike  (1974, 
1975) .  This  limitation  is  due,  in  a  large  part,  to  the  fact  that 
the  structure  of  state  space  models  is  more  general  than  the 
structure  of  time  series  models,  and  the  increase  in  generality 
has  presented  a  significant  challenge  to  the  development  of 
algorithms  that  operate  on  channel  output  data  directly.  Very 
recently,  however,  Van  Overschee  and  De  Moor  (1993)  have  defined  a 
state  space  stochastic  realization  algorithm  which  avoids  the 
computation  of  channel  output  correlation  matrices.  Furthermore, 
this  algorithm  can  be  implemented  using  robust  numerical 
techniques .  The  Van  Overschee-De  Moor  algorithm  was  adopted  in 
Phase  I  to  solve  the  parameter  identification  problem. 

2  .  i  Multichannel _ Detection 

Detection  problems  in  the  context  of  radar  systems  can  be 
postulated  as  hypothesis  testing  problems,  where  a  choice  has  to 
be  made  among  two  or  more  hypotheses.  The  detection  problems 
addressed  in  this  report  involve  the  following  two  hypotheses : 
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H<j:  Desired  signal  is  absent 
:  Desired  signal  is  present 

H0  is  referred  to  as  the  null  hypothesis,  and  H0  is  the  alternative 
hypothesis  .  The  model-based  approach  to  the  multichannel 
detection  problem  is  couched  on  the  assumption  that  the  vector 
random  process  at  the  output  of  the  channels  can  be  represented  as 
the  output  of  a  linear  system  (filter)  under  each  of  the  two 
hypotheses,  and  that  a  unique  parametric  model  corresponds  to  each 
hypothesis.  Furthermore,  the  two  parametric  models  (one  for  each 
of  the  two  hypotheses)  must  be  sufficiently  different  to  allow 
selection  of  the  correct  hypothesis  by  the  evaluation  of  measures 
that  are  sensitive  to  those  differences. 

A  particular  measure  that  has  produced  robust  experimental 
results  in  the  model-based  detection  context  (Metford  and  Haykin, 
1985)  is  the  log-likelihood  ratio  (LLR)  test.  This  test  is  the 
result  of  solving  the  hypothesis  testing  problem  using  the  Neyman- 
Pearson  criterion.  The  LLR  test  in  the  context  of  model-based 
detection  is  calculated  using  the  innovations  sequence  at  the 
output  of  each  of  the  two  linear  filters.  This  presents  practical 
and  implementation  advantages. 

Figure  2-1  illustrates  the  architecture  of  an  on-line 
innovations-based  multichannel  detector.  In  the  case  of  a  radar 
array  system,  each  of  J  radar  receiver  channels  collects  the 
electromagnetic  energy  arriving  at  its  aperture,  and  processes  it 
to  generate  a  discrete-time  random  sequence,  denoted  as  {Xj(n )} , 
which  contains  the  desired  information.  The  J  random  sequences 
(Xj(n)}  are  represented  in  vector  form  as  {^(n)} .  Michels  (1991)  has 

formulated  the  binary  detection  problem  for  multichannel  systems. 
Specifically,  the  null  hypothesis,  H0,  corresponds  to  the  case  of 
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clutter  and  noise  present  in  the  observation  process  {2£(n)},  and  the 
alternative  hypothesis,  H1f  corresponds  to  the  case  of  signal, 
clutter,  and  noise  present  in  the  observation  process  {&(n)} .  That 


is,  the 

detection 

decision  must  be  made  between 

the  following 

models, 

(2-la) 

H0: 

X(n)  «  fi(n)  +  jflt(n) 

n  2:  n0 

(2-lb) 

Hi: 

i(n)  =  a(n)  +  £(n)  +  &(n) 

n  £  n0 

where  n0  denotes  the  initial  observation  time,  (£(n)}  denotes  the 
clutter  process,  {ML(n)}  denotes  all  the  array  channel  noise 
processes,  and  {S(n)}  denotes  the  desired  signal  (target)  process. 
In  the  model-based  approach  pursued  herein,  a  distinct  state 
variable  model  is  associated  with  each  of  the  two  hypotheses,  and 
a  Kalman  filter  is  designed  for  each  model.  Each  Kalman  filter 
processes  the  observation  sequence  {&(n)}  to,  generate  a  vector 
innovations  sequence:  {£(n  |  H0)}  denotes  the  innovations  sequence  at 
the  output  of  the  null  hypothesis  filter,  and  {£(n  I  H 1 )}  denotes  the 
innovations  sequence  at  the  output  of  the  alternative  hypothesis 
filter.  These  innovations  sequences  are  used  in  a  likelihood 
ratio  test  with  a  pre-stored  threshold  to  carry  out  the  detection 
decision . 

As  indicated  in  the  detection  configuration  of  Figure  2-1, 
the  two  filters  can  be  determined  in  real-time  by  processing  the 
observation  sequence  for  a  prescribed  time  interval.  This 
approach  provides  the  most  adaptability,  but  may  present  a  large 
computational  burden  for  some  applications.  It  also  presents 
conceptual  challenges,  such  as  real-time  determination  of  model 
order  for  each  of  the  two  filters.  Alternatively,  the  filter 
design  can  be  carried  out  off-line  for  each  of  the  two  hypotheses, 
and  the  resulting  filter  design  implemented  in  the  real-time 


configuration.  This  alternative  approach  is  less  robust  to 
changes  in  the  operational  environment,  but  requires  a  simpler 
processor  architecture,  which  is  important  in  many  real-time 
applications.  Careful  design  of  the  filters  off-line  using 
adequate  simulated  and  real  data  can  lead  to  acceptable 
performance.  Also,  many  pairs  of  fixed  filters  may  be  designed  to 
cover  distinct  operational  conditions.  The  filter  for  the 
alternative  hypothesis  will  be  of  higher  order  than  the  filter  for 
the  null  hypothesis  because  the  observation  process  for  the 
alternative  hypothesis  has  more  information  (the  signal 
component)  . 


Figure  2-1.  Innovations-based  multichannel  detector  with  on-line 

parameter  identification. 


Michels  (1991)  has  developed  a  likelihood  ratio  calculation 
and  detection  decision  model  which  are  compatible  with  the 
formulation  adopted  herein.  Both  of  these  capabilities  are 


available  at  RL,  and,  where  appropriate,  the  methodology  presented 
in  this  report  is  compatible  with  these  capabilities. 

2 . 2  stattt _ Spacti .  Modal 

The  class  of  multiple-input,  multiple-output  state  variable 
models  can  represent  effectively  the  channel  output  process  for 
radar  applications.  Consider  a  discrete-time,  stationary, 
complex-valued,  zero-mean,  Gaussian  random  process  {^(n)}  defined  as 
the  output  of  the  following  state  space  model  representation  for 
the  system  giving  rise  to  the  observed  process: 


( 2 — 2a ) 

y(n+1)  »  Fy(n)  +  Gu(n) 

n  2  n0 

(2-2b) 

*(n)  =  HHy(n)  +  DHw(n) 

n  £  n0 

(2-2c) 

E&(n0)]  =  Qn 

(2-2d) 

E[x(n0)XH(n0)]  =  P0 

Here  n 

a  n0  denotes  the  initial  time 

(which  can  be  adopted  as 

since  the  system  is  stationary)  .  Also,  ^(n)  is  the  N-dimensional 
state  of  the  system  with  ^(n0)  a  Gaussian  random  vector;  U(n)  is  the 
J-dimensional,  zero-mean,  stationary,  Gaussian,  white  input  noise 
process;  and  ^.(n)  is  the  J-dimensional,  zero-mean,  stationary, 
Gaussian,  white  measurement  noise  process.  The  output  (or 
measurement)  process  (i(n)}  is  also  a  J-dimensional  vector  process. 
Matrix  F  is  the  NxN  system  matrix,  G  is  NxJ  input  noise 
distribution  matrix,  H^1  is  the  JxN  output  distribution  matrix,  D*"1 
is  the  JxJ  output  noise  distribution  matrix,  and  P0  is  the 
correlation  matrix  of  the  initial  state.  All  these  matrices  are 
time-invariant.  Matrix  P0  is  Hermitian  and  positive  definite. 


System  (2-2)  is  assumed  to  be  asymptotically  stable,  which 
means  that  all  the  eigenvalues  of  matrix  F  are  inside  the  unit 
circle.  Also,  system  (2-2)  is  assumed  to  be  reachable  and 
observable,  which  implies  that  the  dimension  N  of  the  state  vector 
(also  the  order  of  the  system)  is  minimal  (Anderson  and  Moore, 
1979)  .  That  is,  there  is  no  system  of  lesser  order  which  has 
identical  input /output  behaviour.  The  output  distribution 
matrices  are  defined  with  the  conjugate  operator  in  order  to  have 
notation  consistent  with  that  of  the  single-output  system  case, 
where  both  H  and  D  become  vectors,  and  nominally  vectors  are 
defined  as  column  vectors. 

The  input  noise  process  correlation  matrix  is  given  as  (all 
matrices  defined  hereafter  have  appropriate  dimensions) 

(2 -3a)  E[u(k)uH(k)]  =  Ruu(O)  =  Q  k  £  n0 

(2 -3b)  E[u(k)uH(k-n)]  »  Ruu(n)  =  [0]  k  £  n0  and  n  *  0 

and  the  output  noise  process  correlation  matrix  is  given  as 

(2 -4 a)  E[^(k)^H(k)]  =  Rww(0)  =  C  k£  n0 

( 2— 4b)  E[w(k)^H(k-n)]  =  Rww(n)  =  [0]  k  >  n0  and  n  *  0 

Notice  that  matrices  Q  and  C  are  Hermitian  (that  is,  =  Q,  and 

CH  =  C)  .  Matrix  Q  .is  at  least  a  positive  semidefinite  matrix  since 
it  is  an  auto-correlation  matrix  (all  the  eigenvalues  of  a 
positive  semidefinite  matrix  are  non-negative) ,  and  matrix  C  is 
assumed  to  be  positive  definite  (this  can  be  relaxed  to  positive 
semi-definite,  but  positive  definiteness  is  more  realistic  since 
in  the  radar  problem  )^(n)  represents  channel  noise  and  other  such 
noise  processes  which  are  independent  from  channel  to  channel) . 


In  the  most  general  form  for  this  model  the  input  and  output 
noise  processes  are  correlated,  with  a  cross-correlation  matrix 

defined  as 

(2 -5a)  E[u(k)^H(k)]  ■  Ruw(O)  -  S  k  2  n0 

(2 -5b)  E[u(k)a£H(k-n)]  «  Ruw(n)  -  [0]  k  2:  n0  and  n  *  0 

In  general,  matrix  S  is  not  Hermit ian.  Both  the  input  and  output 
noise  processes  are  uncorrelated  with  the  present  and  past  values 
of  the  state  process,  and  this  is  expressed  in  terms  of  cross¬ 
correlation  matrices  as 

( ?.  -  6a )  E&(k)uH(k-n)]  =  RyU(n)  =  [0]  k  £  n0  and  n  2  0 

(2-6b)  E|*(k)BH(k-n)] .  Ryw(n)  =  [0]  k  £  n0  and  1120 

The  correlation  matrix  of  the  state  is  defined  as 

(2-7)  E|it(n)it»]  -  Ryyln)  -  P(n)  k  £  n0  and  n  £  0 

It  follows  from  (2-2a)  and  the  above  definitions  that  the  state 
correlation  matrix  satisfies  the  following  recurrence  relation, 

(2-8)  P(n+1)  =  FP(n)FH  +  GQGH  n  s  n0 

In  general,  matrix  P(n)  is  Hermitian  and  positive  definite.  Since 
system  (2-2)  is  stationary  and  asymptotically  stable,  and  since 
matrix  Q  is  positive  definite,  then  the  following  steady-state 
(large  n)  value  exists  for  the  recursion  (2-8) : 


(2-9) 


P(n+1)  =  P(n)  *  P 


Under  steady-state  conditions  Equation  (2-8)  becomes  a  Lyapunov 
equation  for  the  steady-state  correlation  matrix,  denoted  as  P: 

(2-10)  P-FPFh  +  GQGh 

The  conditions  for  steady-state  also  insure  that  the  solution  to 
Equation  (2-10)  exists,  is  unique  (for  the  selected  state  space 
basis),  and  is  positive  definite  (Anderson  and  Moore,  1979). 
Matrix  P  is  unique  for  a  given  state  space  basis.  However,  if  the 
basis  of  the  input  noise  and/or  the  basis  of  the  state  are  changed 
by  a  similarity  and/or  an  input  transformation,  then  a  different 
state  correlation  matrix  results  from  Equation  (2-10) . 

The  correlation  matrix  sequence  of  the  output  process  {&(n)}  is 
defined  as 

(2-iia)  E[x(k)xH(k-n)]  *  Rxx(n)  =  An  V  k  and  n  >  0 

(2-1  lb)  Rxx(-n)  =  Rxx(n)  Vn 

For  a  system  of  the  form  (2-2)  ,  the  correlation  matrix  Rxx(n)  can  be 
factorized  as  follows, 

(2-l2a)  An  =  Rxx(n)  =  n  >  0 

(2-l2b)  An*Rxx(n)  =  rH[Fn-1]HH  =  rH[FH]n‘1H  n<0 

where  F0'1  denotes  F  raised  to  the  (n-l)th  power  and  T  denotes  the 
following  cross-correlation  matrix 


(2-13) 


r  =  E[y(n)*H(n-1 )]  =  Ryx(1 )  =  FP(n)H  +  GSD 


V  n  >  0 


The  correlation  matrix  sequence  factorization  in  Equation  (2-12) 
is  the  key  to  most  correlation-based  stochastic  realization 
algorithms.  The  zero-lag  (n»0)  output  correlation  matrix  is 

(2-14)  Rxx(0)  -  HHP(n)H  +  DhCD  =  Ao 

Matrix  Rxx(0)  is  Hermitian  and  at  least  positive  semidef  inite .  In 

steady-state,  P  replaces  P(n)  in  Equations  (2-13)  and  (2-14) . 

As  can  be  inferred  from  the  above  relations,  the  system 
parameters  {F,  G,  H,  D,  Q,  C,  S,  P.  H  completely  define  the  second-order 

statistics  '(the  correlation  matrix  sequence  {Rxx(n)})  of  the  output 
process,  and  it  is  said  that  system  (2-2)  realizes  the  output 
correlation  matrix  sequence.  Conversely,  the  second-order 
statistics  of  the  output  process  provide  sufficient  information  to 
identify  the  system  parameters,  although  not  uniquely.  Since  the 
output  process  has  zero  mean  and  is  Gaussian-distributed,  the 
second-order  statistics  define  the  process  completely. 

From  the  system  identification  (stochastic  realization)  point 
of  view,  the  problem  addressed  herein  can  be  stated  as  follows: 
given  the  output  data  sequence  {&(n)}  of  system  (2-2)  ,  estimate  a 
set  of  system  parameters  {F,  G,  H,  D,  Q,  C,  S,  P,  T}  which  generates  the 

same  output  correlation  matrix  sequence  as  system  (2-2)  . 
Furthermore,  the  identified  parameter  set  must  correspond  to  a 
system  realization  of  minimal  order  (with  state  vector  y  of  minimal 
dimension) . 

It  is  well  known  (Anderson  and  Moore,  1979)  that  there  can 
exist  an  infinity  of  systems  (2-2)  with  the  same  output 
correlation  matrix  sequence.  The  set  of  all  systems  that  have  the 
same  output  correlation  matrix  sequence  is  an  equivalence  class, 
and  any  two  systems  belonging  to  the  set  are  said  to  be 
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correlation  equivalent  (Candy,  1976) .  For  example,  the  output 
correlation  matrix  sequence  remains  invariant  to  a  similarity 
transformation  applied  to  the  state  vector.  Similarly,  the  output 
correlation  matrix  sequence  remains  invariant  also  to  a  non¬ 
singular  transformation  applied  to  the  input  noise  and/or  to  the 
output  noise.  As  shown  by  Candy  (1976),  the  equivalence  class  of 
correlation  equivalent  systems  is  defined  including  other 
operations  besides  a  change  of  basis . 

Based  on  these  comments,  the  solution  to  the  system 
identification  problem  is  not  unique.  It  is  also  true  that  most 
of  the  possible  system  parameter  solutions  do  not  possess 
desirable  properties.  There  is,  however,  a  solution  which  has 
several  features  of  importance.  This  solution  is  referred  to  as 
the  innovations  representation  for  system  (2-2),  and  is  discussed 
in  Section  2.3. 

In  general,  the  system  matrix  parameters  resulting  from  the 
identification  algorithm  will  be  represented  in  a  different  basis, 
and  should  be  denoted  with  a  different  symbol  (say,  instead  of 

F,  etc.);  nevertheless,  the  same  symbol  will  be  used  in  this 
report  in  order  to  simplify  notation. 

Several  definitions  and  notation  associated  with  the  input 
/output  behaviour  of  system  (2-2)  are  important.  Consider  first 
the  Literm  .(■finite)  controllability  matrix  of  system  (2-2),  C\_i 

this  matrix  is  defined  as  an  NxJL  partitioned  matrix  of  the  form 
(2-15)  »  [g  FQ  FL'1G  ] 

As  is  well-known,  matrix  C\_  has  rank  N  (equal  to  the  system  order) 
for  L  £  N .  The  controllability  matrix  maps  the  input  space  onto 


the  state  space.  Analogously,  the  L-term  observability  matrix  of 
system  (2-2)  is  the  following  JLxN  partitioned  matrix, 


(2-16)  0L 


HH 

hhf 


HhFm 


and  the  rank  of  matrix  <\  is  equal  to  N  also  for  l_£  N.  The 
observability  matrix  maps  the  state  space  onto  the  output  space. 
Classical  realization  theory  for  the  deterministic  case  is  based 
on  the  fact  that  a  block  Hankel  matrix  made  up  of  the  impulse 
response  matrices  (Markov  parameters)  of  a  deterministic  system 
can  be  represented  as  the  product  of  the  observability  and 
controllability  matrices.  That  is, 

(2-17)  HL)U  -  <\C\_ 

where  Hj_l  is  a  JLxJL  deterministic  Hankel  matrix  with  the  impulse 
response  matrix  A(i+j-1)  as  its  (i.j)th  block  element  (a  block  Hankel 
matrix  is  a  matrix  in  which  the  (i,j)th  block  element  is  a 
function  of  i  +  j) .  This  result  follows  from  the  definition  of  the 
impulse  response  matrix  sequence, 

(2-18)  A(n)-HHFn,1G  n*1 


Notice  that  the  factorization  of  the  impulse  response  matrix 
sequence  in  Equation  (2-18)  is  very  similar  to  the  factorization 
of  the  correlation  matrix  sequence  in  Equation  (2-12)  . 

Associated  with  system  (2-2)  is  a  backward  time  model  which 
is  defined  from  the  system  model  (2-2) .  Backward  time  models  play 


a  rola  in  the  formulation  of  a  large  ^lass  of  stochastic 
realisation  algorithm*.  Tha  backward  time  model  for  ayatem  uW) 
is  defined  as  a  diacrata-t ime,  stationary,  complex-valued,  aero- 
mean,  Gaussian  random  procass  with  a  stata  space  ^presentation  of 
tha  form  (Faurre,  1976) 

<2-i9a)  a(n)  -  FHa(n+1)  >^(n) 

(2- 19b)  a(n)  -  I^atn)  ^(n) 


where  a(n)  is  the  N-dimensional  stata  vector,  n )  is  the  N- 
dimensional  input  noise  vector,  and  i£o(n)  is  the  J-dimensional 
output  noise  vector.  Both  noise  vectors  are  uncorrelated  in  time 
(white),  have  mean  equal  to  aero,  and  are  Gaussian-distributed . 
The  L-term  observability  matrix  for  the  backward  system  (2-19)  is 
the  following  JLxN  partitioned  matrix, 


(2-20)  ■ 


r* 

rHFH 


Also  of  interest  is  the  Hermitian  of  'D\_  with  the  block  columns  in 
reversed  order.  That  is, 

(2-21)  (BL  -  *  [FL‘1r  ...  FT  r] 


where  the  dual-point  arrow  over  matrix  indicates  reversal  in 
the  order  of  the  block  columns.  Notice  that  matrix  '3\_  is  like  a 
controllability  matrix  for  the  matrix  pair  (f.  n  in  reverse  block 


oalumn  order.  Thus,  matrix  'J\  \ *  referred  u  herein  xx  the  Lnerm 
reversed  duel  gontraUxbi  \  tty  matrix, 

In  the  oentext  o f  stochastic  realisation  theory,  t  he 
significance  of  the  backward  model  rJlUw*  {from  Equation  (J-'JQ) 
and  the  Henkel  matrix  of  output  correlation  matrices,  aa  shown 
next,  Define  a  ztochaxtio  Henkel  matrix  ,’^A  aa  the  following 

JLxJL  block  matrix, 

Ae  -  >•  \ 

A3  •"  AU1 

\  \  « 

\  *  \ 

\  ^  » 

\*\  "'  *il>\ 

whara  tha  block  elements  (A,)  are  the  elements  of  the  output 
correlation  matrix  sequence,  Equation  (2-12)  ,  It.  follows  from 
Equations  (2-12),  (2-16),  and  (2-22)  that 

(2-23)  9^  -  C\(Cf 

This  equation  is  fundamental  to  stochastic  realisation  algorithms, 
and  allows  the  application  of  classical  deterministic  realization 
algorithms  to  the  stochastic  realization  problem  formulated  with 
output  correlation  matrices.  It  also  provides  insight  into  the 
stochastic  realization  algorithm  presented  in  Section  3.0,  even 
though  the  algorithm  does  not  require  computation  of  the  output 
.‘relation  matrix  sequence. 

Other  important  matrices  in  stochastic  realization  theory 
include  the  JLxJL  "future"  and  "past"  block  correlation  matrices. 
These  matrices  are  the  correlation  matrices  of  future  and  past 
output  block  vectors  defined  as 


(2-22) 


A, 

A. 
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(2-34) 


1?  ■  l(n;jHL-1) 


*(n+i) 


xM-i)  J 


(2-23) 


Xp  •  x(n+l;n+aM) 


»(n*L)  ' 

X(n+L4l) 

\ 

\ 

\ 

.X(n-faL-i) . 


With  these  definitions,  the  future  And  pest  block  correlation 
matrices  are  given  by  the  following  JLxJL  matrices! 


(2-261  SVw  •  ^IWpl 


Ao 

Ac  •• 

■  Am 

<  " 

Aq  « • 

’  AM 

« 

\ 

At-a" 

t  t 

•  A0 

(2-27)  %LX  . 


A0  a,  al-1 


A 


-i  A0  'Va 


Am  Aa>i'‘  Ao  J 


where  and  *:LL  are  the  future  and  past  block  correlation 
matrices,  respectively.  Another  matrix  of  interest  is  the  block 
cross-correlation  matrix  between  the  future  and  the  past,  which  is 
defined  as 


(2-28) 


A^.,  •••  A, 


L  Afl-i  Atl*a  “*  K 


<\A 


Notice  chat  the  block  cross-correlation  Matrix  fyt,Pi  aqua!  to 

the  stochastic  block  Henkel  with  the  block  columns  in  reverse 
order,  aa  indicated  In  Equation  (2-28),  Foe  L  Se  Nr  aquations  ( 2- 
26)-(2~2fl)  define  tha  corralation  structure  of  ayatam  (2-2),  in 
fact,  tha  atochaatic  realisation  algorithm  of  Akaike  (197<),  197*)) 
la  baaad  on  thaaa  block  corralation  matricas. 


2 , 3  fanavitlaBi  Kastif  ntition 


Tha  innovations  rapraaantat ion  ia  a  vary  pawarful  concept  in 
tha  thaory  of  linear  atochaatic  ayatama  dua  to  ita  simplicity  and 
ita  charactarist ica .  Several  texts  and  papers  discuss  this 
concept  in  detail;  in  particular,  Anderson  and  Moore  (1979) 
provide  a  lucid  presentation,  The  discussion  herein  is  adapted 
mostly  from  Anderson  and  Moore  (1979), 

The  innovations  representation  for  a  system  (2-2)  is  a 
discrete-time,  stationary,  complex-valued,  system  of  the  form 

(2-29a)  a(n+1 )  ■  Fa(n)  +  Kf;(n)  nsn0 

(2 -2 9b)  x(n)  -  HHa(n)  +  £(n)  nan0 

(2-2  9c)  alno)-^ 

(2-2 9d)  E[a(n0)aH(n0)]  -  n(n0)  -  n<>  -  [0] 


22 


(2-29*) 

E[a(n)tt>)| .  n(n) 

n  >  ho 

( 2—2  9f ) 

n(n)  ■  n  as  n  ->  ®e 

(2-29g) 

Rxx(n)  *  R**^) 

V  n 

her*  a(n)  la  the  N-dimensional  state,  x(n)  the  J-dimens ional 
output,  and  the  input  process  U(n)}  is  the  innovations  process  for 
•  yatom  (2-2)  ,  That  is,  {fi(n))  is  a  J-dimensional,  sero-mean,  white 
Gaussian  process  with  correlation  matrix  structure  givon  aa 

( 2-30a)  E[c(k)eH(k)]  -  Rmm(O)  •  HHnH  -  \q  •  HHnH  ki  n0 

(2-30b)  E(c(k)£H(k*n)]  •  [0]  k  2  n0  and  ns  0 

The  state  correlation  matrix  H(n)  has  a  steady-state  value  because 

the  system  is  asymptotically  stable  (stationary),  and  the  steady- 
state  value,  n,  is  obtained  as  the  limiting  solution  to  the 

following  recursion 

(2 -3 la)  n(n+1) - FTI(n)FH  +  [Fn(n)H -  H  (Ao  •  HHn(n)H]-1  (Fn(n)H *  qH  n2n0 

<2 -3 ib)  n(n0)  -  a, .  [o] 

Matrix  K  in  Equation  (2-29a)  is  given  as 

(2-32a)  K  ■  [r  -  FHH]  Q'1  -  [r  -  FnH]  (Ao  *  H^H]*1 
(2-32b)  K  ■  GSD  fl'1  ■  QSD  [Ao  •  HHnH]-1 

where  the  second  relation  follows  from  the  definitions  of  T  in 
Equation  (2-13)  and  of  Q  in  Equation  (2-30a)  .  In  the  cases  where 
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tha  inverse  of  the  correlation  matrix  ft  does  not  exist,  its 
pseudoinverse  is  used  instead  in  Equations  (2-31)  and  (.>—32)  . 

Matrices  F,  H ,  A0,  and  T  are  as  defined  for  systtem  (2-2)  . 
That  is,  system  (2-29)  is  related  to  system  (2-2)  ,  In  fact, 
system  (2-29)  as  defined  above  is  the  steady-state  innovations 
representation  for  system  (2-2)  .  This  representation  has  the 
following  important  features. 

(a)  First  and  foremost,  the  correlation  matrix  sequence  of 
(X(n)}  is  .equal  to  the  correlation  matrix  sequence  of 

U(n)},  as  indicated  in  Equation  (2-29g)  .  That  is,  the 
processes  (x(n)}  and  (x(n)}  are  correlation  equivalent. 

This  means  that  the  innovations  representation  is  a 
valid  solution  to  the  system  identification  problem 
defined  herein. 

(b)  Of  all  the  correlation  equivalent  representations  for 
a  given  output  correlation  sequence,  the  innovations 
representation  has  the  smallest  state  correlation 
matrix,  n  (smallest  is  meant  in  the  sense  of  positive 
definiteness;  that  is,  is  smaller  than  112  if  ri2  - 
n,  is  a  positive  definite  matrix) .  This  property  of 
the  innovations  model  is  significant  because  the  state 
correlation  matrix  is  a  measure  of  the  uncertainty  in 
the  state . 

(c)  The  innovations  representation  is  directly  related  to 
the  steady-state  Kalman  filter  (in  the  one-step 
predictor  formulation)  for  system  (2-2) .  In  fact,  the 
steady-state  Kalman  filter  for  system  (2-2)  is 
available  immediately  upon  definition  of  the  steady- 
state  innovations  representation,  and  viceversa. 
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Specif ically,  matrix  K  of  Equations  (2-29a)  and  (2- 
31)  is  the  steady-state  Kalman  gain  of  the  optimal 
one-step  predictor  for  system  (2-2)  .  This  is  true 
provided  that  the  eigenvalues  of  F-KHh  are  stable. 
Thus,  the  innovations  model  is  defined  as  above  for 
all  proqesses  of  the  form  (2-2),  but  the  steady-state 
Kalman  filter  is  defined  only  if  F  -  KHh  is  stable. 

(d)  The  process  (fc(n)}  in  Equations  (2-29)  and  (2-30)  is 
correlation  equivalent,  to  the  innovations  sequence  of 
system  (2-2),  which  is  the  reason  for  referring  to 
system  (2-29)  as  the  "innovations  representation"  for 
system  (2-2)  . 

(e)  The  innovations  model  (2-29)  is  causally  invertible. 
This  means  that  the  present  and  past  of  the  process 
{t(n)}  can  be  constructed  from  the  present  and  past 
values  of  the  output  process  (x(n)} .  The  converse 
statement  is  true  also;  that  is,  any  causally 
invertible  model  is  an  innovations  representation  for 
some  system.  Causal  invert ibility  of  system  (2-29) 
can  be  demonstrated  easily.  From  Equation  (2-29b) , 

(2-33)  E(n)  -  -  HHa(n) -t- x(n) 

Substituting  this  expression  for  £(n)  into  Equation  (2- 
29a)  results  in 

(2-34)  fl(n+1)-[F-KHH]a(n)  +  Kx(n) 


These  relations  demonstrate  the  causal  invertibility 
of  the  innovations  model  (the  input  and  output 
variables  have  traded  places)  . 
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(f)  Matrix  F*KHh  in  the  inverted  innovations  model  is  a 
stable  matrix.  This  follows  from  the  fact  that  the 
matrix  pair  (F,  H)  is  observable,  and  implies  that  the 
Kalman  filter  for  system  (2-2)  is  stable  also. 

(g)  The  transfer  function  of  the  innovations  model  (2-29) 
is  minimum  phase.  This  is  related  to  the  fact  that 
the  innovations  model  is  correlation  equivalent  to 
system  (2-2),  and  second-order  moment  information  (the 
output  correlation  matrix  sequence)  does  not  contain 
any  phase  information. 

(h)  The  innovations  representation  for  a  system  of  the 
form  (2-2)  is  unique.  Given  that  the  innovations 
representation  has  the  same  output  covariance  sequence 
as  system  (2-2) ,  the  fact  that  it  is  unique  eliminates 
searching  for  other  representations  for  system  (2-2) 
with  the  properties  listed  herein. 

(i)  The  innovations  model  (2-29)  can  be  computed  from  the 
output  correlation  matrix  sequence  of  system  (2-2)  . 
This  fact  simplifies  the  parameter  identification 
problem  because  the  set  of  matrix  parameters  that  must 
be  estimated  is  reduced  to  just  five:  (F,  H,  T,  n, 
A0)  (given  these  parameter  matrices,  the  innovations 
covariance,  Q,  and  the  Kalman  gain,  K,  are  obtained 
using  Equations  (2-30a)  and  (2-32a),  respectively). 

All  the  features  listed  above  are  of  relevance  to 
identification  approach  presented  in  Section  3.0  because 
selected  parameter  identification  algorithm  generates 


the 

the 

the 
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innovations  representation  for  the  given  output  correlation  matrix 
sequence,  following  feature  (i) . 

The  backward  model  has  an  associated  backward  innovations 
model  which  is  defined  by  F,  T,  and  the  backward  Kalman  gain. 

Most  of  the  features  (a)-(i)  that  describe  the  forward  innovations 
model  are  valid  also  for  the  backward  innovations  model,  with  a 
notable  exception  of  feature  (b)  ,  which  needs  to  be  replaced  by 
the  following  statement:  For  each  valid  correlation  equivalent 
representation  for  a  given  output  correlation  sequence,  the  state 
correlation  matrix  is  smaller  than  the  inverse  of  the  state 
correlation  matrix  for  the  backward  innovations  model.  More 
specifically,  let  FI^  denote  the  state  correlation  matrix  for  the 
backward  innovations  model  in  steady-state  conditions,  and  let  £ 
denote  the  state  correlation  matrix  for  any  valid  correlation 
equivalent  representation  of  an  output  correlation  sequence. 

Then,  n;1  -  £  is  a  positive  definite  matrix.  This  result  provides 

an  upper  bound  for  the  state  correlation  matrix  of  a  correlation 
equivalent  representation,  and  can  be  combined  with  the  lower 
bound  established  by  property  (b)  of  the  forward  innovations  model 
to  give 

(2-35)  n  £  £  £  n^1 

As  before,  the  inequality  between  two  matrices  is  intended  in  the 
sense  of  positive  semi-definiteness  of  the  matrix  difference. 
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3 . 0  MULTICHANNEL  SYSTSM  IDENTIFICATION 


Identification  of  the  model  parameter  matrices  (F,  H,  T,  IT,  A<>}  is 
carried  out  using  the  algorithm  of  Van  Overschee  and  De  Moor 
(1993),  extended  to  the  case  of  complex-valued  data.  The  Van 
Overschee-De  Moor  algorithm  is  based  on  the  predictor  space 
concept  of  Akaike  (1974;  1975),  the  correlation  equivalence 
results  obtained  by  Faurre  (1976),  and  the  balanced  stochastic 
realization  approach  of  Arun  and  Rung  (1990).  The  algorithm 
approach  is  presented  herein  from  a  viewpoint  which  is  different 
from,  and  simpler  than,  the  presentation  given  by  Van  Overschee 
and  De  Moor  (1993)  . 

3 . 1  Output  Data-Baaed  Algorithm 

In  comparison  with  alternative  stochastic  realization 
techniques,  the  Van  Overschee-De  Moor  algorithm  adopted  herein  has 
several  advantages  for  multichannel  detection  applications,  as 
listed  next. 

•  Reduced  dynamic  range  with  respect  to  algorithms  which 

require  generation  of  the  output  correlation  matrix 
sequence  (correlation  matrices  are  estimated  as  sums  of 
products  of  the  data  sequence  elements,  which  increases 
the  dynamic  range)  .  As  such,  the  algorithm  can  be 

viewed  as  a  "square-root"  algorithm. 

•  Identifies  the  parameters  for  a  model  in  the  state-space 
class,  which  is  more  general  than  the  time  series  class. 

•  Belongs  to  a  class  of  algorithms  referred  to  as 
"subspace  methods."  Subspace  methods  involve  the 
decomposition  of  the  space  spanned  by  the  output  process 
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into  two  orthogonal  subspaces:  one  subspace  is  the  space 
spanned  by  the  "desired  component,"  and  the  other 
subspace  is  spanned  by  the  "noise  component."  The  MUSIC 
algorithm  (Schmidt,  1979;  1981),  for  example,  also 

belongs  to  the  class  of  subspace  methods. 

•  An  approximately  balanced  (in  the  stochastic  sense) 
state  space  realization  is  generated,  thus  providing  a 
built-in  and  robust  mechanism  for  model  order  selection. 

•  Identifies  the  innovations  representation  of  the  system, 
and  generates  the  Kalman  gain  directly,  without  having 
to  solve  a  nonlinear  discrete  matrix  Riccati  equation. 

•  Approach  differs  from  others  in  that  the  states  of  a 
Kalman  filter  for  the  given  sequence  are  identified 
first,  and  then  the  model  parameters  are  estimated  via 
least-squares . 

•  Implementation  of  the  algorithm  involves  the  QR 
decomposition  and  the  quotient  SVD  (QSVD;  also  known  as 
the  generalized  SVD),  which  are  stable  numerical 
methods.  Furthermore,  the  QSVD  is  applied  to  matrices 
of  small  dimensions. 

An  algorithm  for  implementing  the  QSVD  is  given  in  Appendix  B  for 
the  specific  conditions  presented  in  this  section. 

Consider  the  channel  output  sequence  {i(n)} .  For  simplicity, 
let  the  initial  time  no  =  0.  This  can  be  done  without  loss  of 

generality  because  the  system  is  stationary.  Now  define  a  block 
Hankel  matrix  X0  ^  with  output  sequence  vectors  assigned  as  block 

elements  according  to  the  rule  X0  |_.i(i.j)  =  &(i+j-2);  that  is, 
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'  m 

am 

X(2)  ••• 

S(2) 

&(3)  •  •  • 

*(M) 

(3-1) 

*0,1-1  » 

*(2) 

i(3) 

X(4)  ••• 

*(M+1) 

a(L)  x(L+1)  •  •  • 

*(L+M-2) 

Here  the  first  subscript  denotes  the  time  index  of  the  first 
element  of  the  first  row,  and  the  second  subscript  denotes  the 
time  index  of  the  first  element  of  the  last  row.  Matrix  Xq  |_.i  has 

JL  rows  and  M  columns,  with  M  » JL,  and  JL  >  N  (recall  that  N  is 
the  system  order  and  J  is  the  number  of  channels) .  The  block  row 
dimension,  L,  must  be  selected  so  that  L  N  +  1  .  Similarly,  define 
another  JLxM  block  Hankel  matrix  Xl2l-i  with  output  sequence  vectors 
assigned  as  block  elements  according  to  the  rule  X^l-iO'D  =  X(i+j*2+L); 

that  is, 


- 

m  *(L+i) 

*(L+2) 

...  x(L+M-1) 

S<L+1)  i(L+2) 

X(L+3) 

•  •  •  x(L+M) 

(3-2) 

*L.2L-1  " 

S(U2)  js(L+3) 

*(L+4) 

. ..  &(L+M+1) 

)  s(2L) 

X(2L+1) 

...  *(2L+M-2) 

Matrices  Xq  and  X|_  2L-1  rePresen^  the  "past"  and  the  "future", 
respectively,  of  the  output  process.  Also  let  X~  denote  the 
vector  space  spanned  by  the  past  of  the  process  (X.(n)},  and  X + 
denote  the  vector  space  spanned  by  the  future  of  the  process . 
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The  algorithm  is  based  on  the  decomposition  of  the  process 
future  (as  represented  by  matrix  XL2L..j)  into  two  orthogonal 

subspaces  (herein  orthogonality  is  invoked  under  the  standard 
unitary  inner  product  for  complex  vector  spaces,  with  identity 

metric)  .  In  such  a  decomposition,  one  subspace  is  the  space 
spanned  by  the  process  past,  X~  (as  represented  by  matrix  X0l.-|)> 

and  the  second  subspace  is  the  space  spanned  by  the  noise  process. 
Let  *W  denote  the  space  spanned  by  the  noise  process  {^(n)} .  Then, 
the  desired  decomposition  of  X*  is  as  follows: 

(3-3)  X*  -  XT  ©  •W 

where  ©  denotes  the  direct  sum,  and  X~  1  *W  (since  the  present 

state  and  the  present  measurement  noise  are  uncorrelated) .  In 
matrix  notation,  the  desired  decomposition  of  X|_ 2l-i  exPressed  as 

(3-4)  ^L,2L-1  =  *FiiP  +  ^FIP 

where  the  JLxM  matrix  Xp(lP  is  the  projection  of  the  row  space  of 
X|_,2L-1  future)  onto  the  row  space  of  X0L..|  (the  past),  and  the 

JLxM  matrix  XFlP  is  the  projection  of  the  row  space  of  X(_2L.1  onto 
the  complement  of  the  row  space  of  X0  ^ .  Akaike  (1974  ;  1975)  has 
demonstrated  that  since  the  order  of  the  state  space  model  is  N, 
the  projection  of  the  future  onto  the  past’  is  an  N-dimensional 
subspace  of  the  M-dimensional  space  to  which  the  rows  of  X^i-i 
belong.  Thus,  the  rank  of  matrix  XF||P  is  equal  to  the  dimension  of 

the  subspace  spanned  by  the  projection  of  the  future  onto  the 
past.  Furthermore,  the  structure  of  this  subspace  (and  of  its 

matrix  representation)  determines  the  characteristics  of  the  state 
space  model  (such  as  model  order)  .  Analogously,  matrix  Xpip 

determines  the  characteristics  of  the  noise  subspace. 


The  decomposition  (3-4)  can  be  carried  out  using  a  matrix 
operator  referred  to  as  a  projector  (Pease,  1965) .  Let  and  S2 
denote  two  orthogonal  subspaces  of  5  such  that  «5i  ©  S2  =  *5*  A 
projector  of  5  onto  is  a  matrix  P1  such  that 

(3-5a)  Pi¥i=¥i  v¥i€^i 

(3-5b)  P^-Q  Vy2eJ2 


Projectors  can  be  defined  also  as  operating  on  row  vectors, 
instead  of  on  column  vectors.  The  property  which  characterizes 
projectors  is  idempotency  (that  is,  P  is  a  projector  if  and  only 
if  P2  =  P)  . 

Let  V  denote  the  M-dimensional  subspace  defined  by  the  L  rows 
of  XL2L..,  (recall  that  L«M),  and  let  P.  denote  the  MxM  projector 
of  V  onto  the. subspace  X~ .  It  follows  that 


(3-6) 


XL2L-lP- 


>nP 


Thus,  availability  of  the  projector  P.  allows  the  decomposition  of 
the  future  data  matrix  because  XFj_p  can  be  determined  from 

Equations  (3-4)  and  (3-6) .  Projector  P_  is  determined  from  matrix 

Xq.L-1  as 


(3-7) 


*0,L*1  *0,1*1 


*0.1*1  (*0,L-1*0,L*l)  *0.L-1 


This  projector  will  decompose  the  future  data  matrix  into  the 
desired  components.  However,  Equation  (3-7)  imposes  a  large 
computational  burden,  and  furthermore,  it  effectively  involves  the 
calculation  of  the  output  correlation  matrix  sequence  and  of  the 
inverse  of  a  large  matrix  with  correlation  matrix  sequence 
elements  as  its  block,  elements.  Fortunately,  the  QR  decomposition 
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can  ba  applied  to  determine  the  subspace  decompoait ion  and 

also  to  determine  the  projector  P„,  if  required.  The  OR 

decomposition  is  a  computationally  efficient  and  numerically 
robust  approach  to  address  this  problem  (Dongarra  et  al,,  1979). 


Consider  now  the  block  Hankel  data  matrix  X0t2L>l  •  This  matrix 
is  a  2JLxM  block  column  matrix  made  up  of  a  concatenation  of  the 
past  and  future  Hankel  matrices , 


(3-8) 


X 


0,21-1 


*0,1-1 

.  XL,2L-1  - 


Now  apply  the  Hermitian  operator  to  a  "normalized"  matrix  X0t2L-l ' 
and  carry  out  a  QR  decomposition  on  this  matrix  (the  normalization 
factor  is  required  to  avoid  increase  in  dynamic  range  and  to  match 
the  formulation  which  is  based  on  the  correlation  matrix 
sequence) .  This  results  in 


(3-9a) 
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Matrix  Q  is  an  MxM  unitary  matrix,  submatrices  QA  and  QB  are 
dimensioned  MxJL,  and  submatrix  Qq  is  dimensioned  Mx(M-2JL), 

Matrix  R  is  a  2JLx2JL  upper-triangular  matrix  with  rank  equal  to 
the  rank  of  matrix  X0i2L-1  •  All  t^ie  submatrices  of  R  are 
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y  y 

dimensioned  JLxJL,  and  submatrica*  RA  and  Rq  are  also  upper- 
triangular.  Since  matrix  Q  ia  unitary,  tha  fallowing  relations 
ara  truat 


(3-10)  QQh  -  QaQJ  +  Q8Qj  +  Q0Qq  •  lM 
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QaQb 
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Conaidar  now  tha  conjugate  transpose  of  Equation  (3-9)  ,  after 
eliminating  Qq  since  it  is  multiplied  by  a  zero-valued  matrix/ 

that  is, 


(3-12) 


X0_,2L-1  ^  1 

X0.L-1 

■  *A 

(0)  ' 
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-  Xl,2L-1  . 
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The  following  two  equations  are  obtained  immediately  from  the 
partitioning  in  Equation  (3-12), 

(3-13)  ^  -  RAQj 

(3-h)  ^  -  RbQa  ♦  RcQg 

Equation  (3-13)  is  a  QR  decomposition  of  X0ti_.i  (recall  that  RA  is 

lower  triangular),  and  Equation  (3-14)  is  a  subspace  decomposition 
°f  X|_2L.i  •  As  shown  next,  (3-14)  is  the  desired  decomposition  of 
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XL(|L%1 ,  Tht  projector  P.  ia  determined  from  Equations  (3-7)  and  (3- 
13)  at 

(3-15)  P_  -  QAQj 

Zt  alao  follows  that  tha  projactor  of  V  onto  W  (tha  noisa 
•ubapaca)  ia 

(3-16)  P*>  ■  QgQgj 

With  Equations  (3-11)  and  (3-15)  it  ia  easy  to  demonstrate  that 

u-17)  \aL.,p.  .  VTT  (  r.q” ]  .  xFl# 

Similarly,  it  follows  from  Equations  (3-11)  and  (3-16)  that 
(3-18)  XL,2L-lP^  "  [  RcQbI  *  ^FiP 

This  demonstrates  that  Equation  (3-14)  is  the  decomposition  of  the 
future  onto  the  past  and  onto  the  noise  orthogonal  subspaces. 

The  information  of  the  projection  of  the  future  onto  the  past 
is  contained  in  matrix  Rg.  Specifically,  the  rank  of  Rg  is  equal 

to  the  order  of  the  state  space  model  representation  for  the 
future-to-past  interface,  and  the  column  space  of  Rg  is  equal  to 

the  column  space  of  the  observability  matrix  for  the  state  space 
model  (Van  Overschee  and  De  Moor,  1993) . 

At  this  point  in  the  development  it  is  convenient  to  continue 
the  decomposition  of  the  Q  and  R  matrices  in  Equation  (3-12)  in 
order  to  isolate  as  much  as  possible  the  structure  of  the 
orthogonal  subspaces.  To  that  end,  consider  Equation  (3-12),  and 
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carry  out  a  further  partitioning  of  tha  QR  decomposition  matrices 
aa  follows  (tha  dimensions  of  tha  matricaa  on  tha  right-hand-side 
of  Equation  (3-19a)  are  given  in  (3-19b)  and  (3-19c)); 


(3-19a) 
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From  Equations  (3-12)  and  (3-19)  it  follows  that  the  JLxJL  matrices 
Rb  and  Rq  are  defined  with  the  following  partitions: 

(3-20)  Rg  * 


R33  [0] 

R43  R44  - 


(3-21)  Rc  = 


*31  *32 

.  R>11  Ra9  . 
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Rtftr  to  th*  partitioning  in  Equation  (3-19)  and  define  two  other 
partitioned  matrices  as 

(3-22)  Rq  ■  [  R41  R^j  R43  ] 

(3—23)  Rg  ■  R44 

Matrix  R0  is  J(L*1)xJ(L+1 ),  and  matrix  Rg  is  J(L-1)xJ(L-1) .  Now  carry 

out  two  QSVDs  •-'n  these  matrix  pairs  as  detailed  in  Appendix  B, 

One  QSVD  is  app  Led  to  the  matrix  pair  Rg  and  Rc  to  obtain 

(3-24)  Rg  -  UlSlYl 

(3-25)  Rc  -  VlTlYl 

The  second  QSVD  is  carried  out  on  the  matrix  pair  RD  and  Rg,  and 
results  in 

(3-26)  Rq  -  UL.,SL.,Yg., 

(3-27)  Rg  .  VL.,TL.tYf.t 

In  these  two  QSVDs,  matrices  Ul-i,  Ul,  V^,  and  V[_  are  unitary,  and 

matrices  Y|__i  and  YL  are  square  and  non-singular  (Appendix  B)  . 

Also,  the  subscripts  (L  or  L-1 )  correspond  to  the  term  index  of  an 
associated  observability  matrix  defined  as  in  Equation  (2-16). 
That  is,  the  following  two  results  are  true: 

•  The  column  space  of  matrix  RD  is  the  same  as  the  column 
space  of  0|_.i  .  This  result  follows  from  two  facts: 
first,  for  L-1  >  N  the  observability  matrix  maps  the 
state  space  onto  the  output  space;  and  second,  the 
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dtcompoa it  ion  on  Equation  (3-19)  indicate?  a 

decomposition  of  the  form  (3-14)  for  the  block  Hankal 
matrix  t  which  consists  of  tha  last  L*1  block 

rows  of  Xqjl.i  in  Equation  (3-19)  . 

•  Tha  column  space  of  matrix  Rg  is  tha  same  as  tha  column 
space  of  0|,.  As  in  tha  preceding  argument,  this 

follows  from  tha  mapping  property  of  the  observability 
matrix  and  from  the  decomposition  of  matrix  X^l-i  in 

Equation  (3-14)  . 

Consider  the  cases  where  the  matrix  pairs  (Rg,  Rq)  and  (Rq,  Rg) 
form  concatenated  matrices  of  full  rank  (see  Appendix  B)  ,  which 
are  the  most  likely  cases  in  practical  situations  involving  random 
data.  In  those  cases  matrix  Sl.i  is  rectangular  _h  2J  more  rows 

than  columns,  and  is  zero  except  possibly  along  the  main  diagonal. 
The  elements  along  the  main  diagonal  of  Si_.i  are  real-valued,  with 

value  bound  between  unity  and  zero,  and  arranged  in  order  of 
decreasing  magnitude.  Matrices  S|_,  T|_.i,  and  T|_  are  square  and 
diagonal.  The  diagonal  elements  of  S|_  are  real-valued,  with  value 

bound  between  unity  and  zero,  and  arranged  in  order  of  decreasing 
magnitude  also.  The  diagonal  elements  of  both  and  Tl  are  also 

real-valued  and  with  value  bound  between  unity  and  zero.  However, 
the  diagonal  elements  of  these  two  matrices  are  arranged  in  order 
of  decreasing  magnitude.  In  pairs,  the  diagonal  elements  of 
and  Tl.i  are  referred  to  as  singular  value  pairs  of  matrices  Rq  and 
Rg.  Likewise,  the  diagonal  elements  of  S|_  and  T[_  are  the  singular 
value  pairs  of  matrices  Rg  and  Rq. 

The  value  of  the  diagonal  elements  of  matrices  and  S|_  are 

indicative  of  model  order.  In  fact,  when  the  data  is  the  output 
of  a  system  of  order  N,  only  the  first  N  diagonal  entries  are  non¬ 
zero  in  matrices  and  S|_  (model  order  determination  is 
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discussed  further  in  Section  3.2).  As  discussed  in  Appendix  B, 
for  every  zero-valued  diagonal  entry  in  matrix  Sl  there  is  a 
corresponding  unity-valued  diagonal  element  in  matrix  T^.  The  same 
relationship  is  true  for  matrices  S^-i  end  .  Thus,  for  an  N-th 
order  model  the  two  pairs  of  S(,j  and  T(,j  matrices  have  a  natural 
partition  along  the  main  diagonal  corresponding  to  the  first  N 
entries.  Specifically, 

[0] 

*JL-N  - 

[0] 

^JUJ-N.JL-J-N  - 
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Now  define  block  column  partitions  in  matrices  U(.),  V(.),  and  Y(.j  to 


correspond  with 

results  in 

the 

partitions  in  Equations 

(3-28) - (3-31) . 

This 

(3-32)  UL  =  [ 

yO) 

u«]  uL,=[u<:,' 

i 

(3-33)  VL  =  [ 

V(D 

v"]  vL,.[v<;» 

vg] 

(3-34)  Yl  =  [ 

y(1) 

tl 

Y«]  yl, 
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(3-34) 


All  submatrices  with  superscript  1  have  N  columns.  Partitioning 
of  the  S(.),  T(,)f  U(.),  V(.),  and  V(.)  matrices  does  not  provide  a 

theoretical  advantage  or  an  enhanced  insight;  it  does,  however, 
reduce  the  computational  burden. 


The  two  QSVDs  were  introduced  to  extract  the  structure  and 
subspace  information  available  in  the  R  matrix  (and  submatrices) 
of  the  QR  decomposition.  Substitution  of  the  QSVD  results  into 
the  corresponding  partitions  in  Equations  (3-12)  and  (3-19)  allows 
appreciation  of  this  structure.  Carrying  this  out  leads  to  the 
following  expressions: 


(3-35) 
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(3-36) 
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Both  of  these  equations  exhibit  a  subspace  decomposition  of  the 
respective  block  Hankel  data  matrix,  and  in  each  equation  the 
information  of  the  structure  of  the  two  orthogonal  subspaces  is 
contained  in  the  partitioned  matrix  involving  the  S(.)  and 

matrices.  In  Equation  (3-36)  the  partitions  in  the  matrices  are 
not  emphasized  because  the  dimensions  of  the  individual  partitions 
are  not  compatible  (as  they  are  in  Equation  (3-35)) .  Of  course, 
overall  matrix  dimensions  are  compatible. 


Given  the  decompositions  in  Equations  (3-35)  and  (3-36),  it 
remains  to  develop  the  procedure  that  relates  these  decompositions 
to  the  innovations  model  parameter  matrices.  This  is  done  using 
orthogonal  projections  in  random  vector  spaces. 
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Consider  the  orthogonal  projection  of  the  future  block  vector 
Xp  onto  the  past  block  vector  Xp  (recall  Equations  (2-24)  and  (2- 

25) ) .  This  projection  is  the  minimum  variance  estimate  of  Xp  given 
Xpr  which  is  also  the  conditional  mean  of  Xp  given  Xp.  That  is, 

(3-37)  Xp  -  EtXpIXp]  -  EtxpXpHEtaXp])'1  Xp 


Using  Equations  (2-26)  and  (2-28)  this  can  be  expressed  compactly 
in  terms  of  the  observability  and  reversed  dual  controllability 
matrices  as 

(3-38)  Xp  =  ^,L^:L,L*P  =  ^:L,  L*P 


Suppose  a  minimum  variance  estimate  is  sought  for  each  one  of  the 
columns  of  the  data  block  Hankel  matrix  X|_j2l-i ^  which  represents  the 

future.  Then,  concatenating  M  such  estimate  equations  into  a 
single  matrix  estimate  equation  leads  to 

(3-39)  X|_,2l-1  =  ^p.L, L^0,L-1  =  OJ-L 


The  NxM  matrix  ZL  is  very  important,  and  deserves  to  be  defined 
directly,  as  in  Equation  (3-40)  next, 

(3-40)  ZL  =  ■^P:L,  L^0,L-1 

Equation  (3-39)  states  that  the  minimum  variance  estimate  of  the 
output  (the  columns  of  matrix  is  a  linear  function  of  the 

columns  of  matrix  ZL.  Recall  that  the  observability  matrix  maps 

the  state  vector  into  the  output,  and  that  the  states  of  a  Kalman 

filter  are  minimum  variance  estimates  of  the  states  of  the  linear 
system  for  which  the  filter  is  designed.  Thus,  the  columns  of  ZL 


are  states  of  a  Kalman  filter  for  the  system  to  be  identified. 
This  result  is  instrumental  to  the  algorithm. 

Equations  (3-39)  and  (3-40)  involve  correlation  matrices, 
which  is  undesired  because  of  the  computational  burden  associated 
with  their  calculation  and  also  because  of  the  increase  in 
numerical  precision  (dynamic  range)  associated  with  computations 
involving  correlation  matrices.  It  is  possible  to  convert  these 
equations  to  "square  root"  form  by  substituting  estimates  of  the 
correlation  matrices  calculated  using  channel  output  data.  For 
sufficiently  large  values  of  M  the  reversed  stochastic  Hankel 
matrix  and  the  past  block  correlation  matrix  are  approximated 
effectively  by  the  biased  correlation  matrix  estimators  using  the 
channel  output  data,  and  it  is  simple  to  demonstrate  that  such 
estimates  can  be  represented  in  terms  of  the  output  data  block 
Hankel  matrices.  Specifically,  for'  large  M, 

(3-41)  ^P:L, L  *  ■jjj  k).L-lX0,L-1  *  RARA 

(3-42)  \iL  .  Jj.  X^Xl-l  -  wW!  -  <VK 

where  Equations  (3-13)  and  (3-35)  have  been  applied.  Equation  (3- 
42)  suggests  that  the  reversed  stochastic  Hankel  matrix  can  be 
factorized,  as  in  deterministic  realization  problems  (Zeiger  and 
McEwen,  1974),  to  obtain  the  observability  and  the  reversed  dual 
controllability  matrices, 

(3-43)  q.  -  YlSJ*  -  y(l1)(s[1))1/z 

0-44)  \  -  s^u^rJ  -  (s^)1/2(u^)hrX 
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Now  the  matrix  of  Kalman  filter  states,  ZL,  can  be  determined  as 
(using  Equations  (3-13),  (3-40),  (3-41),  and  (3-44)), 


(3-45,  Z,  =  S’XdJ  =  (8‘T H")HQ»  -  (sfriufT 


H 


Q1 


q! 


H 


This  is  a  key  result  of  the  algorithm. 


In  order  to  calculate  the  system  parameter  matrices  F,  T,  and 
H  it  is  necessary  to  determine  two  additional  filter  state 
matrices:  Z l+1  and  Wj_.  Matrix  ZL+1  can  be  thought  of  as  a  "shifted" 
version  of  matrix  ZL;  in  fact,  the  columns  of  matrix  ZL+1  are 
obtained  from  the  .columns  of  matrix  ZL  via  a  Kalman  filter  (or  an 
innovations  model,  Equation  (2-34)).  Determination  of  matrix  ZL+1 
requires  steps  identical  to  those  in  the  derivation  of  matrix  ZL. 

The  main  difference  is  that  the  matrices  involved  correspond  to 
the  (L+1)-terrri  observability  and  reversed  dual  controllability 

matrices.  And  consequently,  the  results  of  the  QSVD  of  matrices 
Rq  and  Rg  are  utilized.  Repeating  the  steps  in  the  derivation  of 

matrix  ZL  leads  to  the  following  result: 


(3-46) 


(s<,,)','2(xn>),Ym  s'’»(u‘’»)H 


Q1 

Q! 

Q1 


where  the  underbar  denotes  that  matrix  is  obtained  from  matrix 
by  deleting  the  last  block  row. 

If  the  procedure  to  determine  matrix  Zj_  is  followed  departing 
from  the  orthogonal  projection  of  the  past  block  vector  onto  the 
future  block  vector  then  the  resulting  matrix,  denoted  as  W(_, 

is  a  matrix  of  Kalman  filter  states  for  the  backward  Kalman  filter 
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(the  Kalman  filter  for  the  backward  system)  .  The  form  of  this 
backward  filter  states  matrix  is 


,3-47)  WL  -  (S<”H  S<»K”f  Tf>(vf>r  ] 


The  main  difference  between  this  expression  and  the  expression  for 
matrix  Equation  (3-45),  is  that  Equation  (3-47)  includes 

matrices  TL  and  VL  directly. 

Given  the  Kalman  states  matrices  ZL,  ZL+J,  and  WL,  the  system 
parameter  matrices  F,  T,  and  H  can  be  determined  as  least-squares 
solutions  to  linear  systems  of  equations  in  noise.  This  is  a 
result  of  the  relationship  with  these  matrices  and  the  forward  and 
backward  Kalman  filter  for  the  channel  output  sequence  (£(n)} .  The 
procedure  is  described  next . 

Since  the  columns  of  the  matrices  ZL  and  ZL+1  are  Kalman 
filter  states,  it  follows  that 

(3-48)  Zu,  «  FZt  +  Zy 

where  Zy  is  a  matrix  of  residuals  orthogonal  to  (Kalman  filter 

residuals  are  orthogonal  to  the  state  estimates) .  A  least-squares 
estimate  of  F  is  obtained  from  Equations  (3-45),  (3-46),  and  (3- 

48)  as 

(3-49.)  F  . 
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(3-49b) 


(st,>r/2w’,),Y«:»s«’.>(u<:.'1)Hu<”(s[,»r 
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Alternative  formulas  for  matrix  F  can  be  defined  based  on  several 
relationships  (observability;  controllability;  backward  model; 
etc.)  that  are  valid  for  the  system  matrix.  Each  distinct  formula 
presents  different  numerical  precision  and  computational 
requirements .  This  issue  defines  an  important  set  of  trade-offs 
for  investigation  in  Phase  II. 


The  output  equation  for  the  Kalman  filter  leads  to  the 
following  matrix  relation, 

(3-50)  Xl-H^Zl  +  Z* 

where  matrix  Xl  is  the  first  block  row  of  matrix  XL2L-  and  Zyy  is  a 
matrix  of  residuals  orthogonal  to  Zj_.  From  the  definition  of  Xl 
and  Equations  (3-12)  and  (3-19), 


(3-511  ^-(sd)  *(L+1)  •••  »(L*M-1)]  -  [  Rj,  R*  R*  (01] 


From  Equations  (3-45),  (3-50),  and  (3-51)  the  least -squares 

estimate  of  HH  is 


(3-52»>  Hh  ■  x,_z[  -  XLZ^( z,_z^ )'’ 


HH  -  [  R31  fl*K’(s<”)' 


1/2 
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(3-52b) 


This  is  a  simple  expression  and  it  involves  matrices  of  relatively 
low  dimensionality.  An  alternative  derivation  for  Equation  (3- 
52b)  departs  from  Equation  (3-43)  and  capitalizes  on  the  fact  that 
matrix  H  occupies  the  first  J  rows  of  the  observability  matrix. 


A  least-squares  estimate  for  T  is  obtained  in  a  manner 

U 

analogous  to  the  solution  for  H  obtained  above.  This  is  based  on 
the  fact  that  is  the  output  measurement  matrix  for  the  backward 

system.  Thus,  the  output  equation  for  the  backward  Kalman  filter 
leads  to  the  following  matrix  relation, 

(3-53)  Xl.,  -  ^Wl  +  Zv 

where  matrix  is  the  last  block  row  of  matrix  X0L..,,  and  Zv  is  a 
matrix  of  residuals  orthogonal  to  WL.  From  the  definition  of  X^ 
and  Equations  (3-12)  and  (3-19), 


(3-54)  X^  -  [  *d*1)  X(L)  •••  i(L+M-2)  ]  -  [  R2i  R22  t°l 


Then,  based  on  Equations  (3-47),  (3-53),  and  (3-54),  the  least- 

squares  estimate  of  r*  is  obtained  as 


(3-55a)  Th  -  XL.,W*  -  X^WPKwr)'1 
(3-55b)  rH.[R2,  R2]u[1)(s(M),/2 


This  expression  is  analogous  to  Equation  (3-52)  .  Just  as  in  the 
case  for  H  ,  an  alternative  derivation  of  Equation  (3-55b)  is 
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possible  based  on  the  fact  that  matrix  T  occupies  the  last  J 
columns  of  the  reversed  dual  controllability  matrix,  Equation  (3- 
44)  . 


Alternative  formulas  can  be  defined  also  for  matrices  H  and  T 
based  on  the  various  system  relationships  that  involve  these 
matrices.  This  also  constitutes  an  important  set  of  trade-offs 
for  Phase  II. 

Notice  that  the  Q(.)  matrices  do  not  appear  in  the  final 

expressions  for  the  matrix  parameters.  The  QR  decomposition  is 
fundamental  to  the  algorithm,  but  only  the  R(.)  matrices  have  to  be 

calculated  and  stored.  This  is  a  very  important  feature  of  the 
algorithm  because  one  dimension  of  Q(,)  is  very  large  (M),  and  the 

manipulation  of  these  matrices  would  involve  significant  storage 
and  computational  requirements. 

Determination  of  the  remaining  matrix  parameters  for  the 
innovations  model  (2-29)  is  described  next.  Consider  first  the 
steady-state  .correlation  matrix  of  the  innovations  model  state,  n. 

This  correlation  matrix  is  equal  to  the  correlation  matrix  of  the 
Kalman  filter  state  (Anderson  and  Moore,  1979)  .  Therefore,  a 
robust  estimator  for  FI  is  based  on  the  columns  of  matrix  Z l, 

(3-56)  n  =  ZjZl  =  s(L1} 


It  turns  out  that  the  backward  filter  states  also  lead  to  the  same 
result, 

0-57)  =  wLw^  = 

A  system  model  such  that  the  forward  and  backward  correlation 
matrices  are  both  diagonal  and  equal  is  said  to  be  in  balanrpd 


poardinitti  (in  th«  stochastic  sonsa) .  Balancsd  coordinatts  allow 
•ffactiva  modal  ordar  salactiort  and/or  raduction. 

Tha  zaro-lag  output  corralation  matrix  is  obtainad  diractly 
from  tha  output  saquanca  as 

(3-58)  Ao  ■  -J-  £ 

Ny 

(3-59)  Nj  ■  M  +  2L-1 

hara  Nj  is  tha  total  number  of  output  data  vectors  (length  of  the 

output  sequence)  used  in  the  algorithm.  The  innovations 
correlation  matrix  is  obtained  from  Equation  (2-30a), 

(3-60)  Cl-Ao*HHnH 

Finally,  the  one-step  prediction  filter  (Kalman)  gain  is  obtained 
from  Equation  (2-32a)  as 

(3-61)  K  ■  [r  -  FTIH]  il'1  m  [r  -  FT1H]  [Ao  -  HHnH]‘1 

which  completes  the  model  parameter  identification  algorithm. 

3 . 2  Modal  Order  Determination 

Model  order  determination  is  a  necessary  decision  for  any 
identification  algorithm  in  applications  where  the  true  order  of 
the  system  generating  the  channel  output  data  is  unknown,  or  where 
the  true  process  generating  the  data  may  not  belong  to  the  model 
class  adopted  to  represent  the  data.  In  the  second  case  the  model 
generated  by  the  algorithm  is  a  "representation  model,"  as  opposed 
to  a  "physical  model"  (a  model  based  on  accurate  analyses  of  the 
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underlying  physical  processes) .  Determination  of  the  model  order 
is  always  a  difficult  problem,  and  the  solution  is  rarely  clear- 
cut.  The  Van  Overschee-De  Moor  identification  algorithm  does  have 
several  strong  features  that  lead  to  robust  model  order 
estimation.  Principally,  the  algorithm  identifies  the  model 
parameters  of  the  innovations  representation  for  the  multichannel 
process  in  balanced  coordinates.  Model  order  determination  is  an 
important  issue  and  deserves  detailed  analysis  during  Phase  II. 

The  prime  mechanism  for  model  order  selection  in  the 
algorithm  is  examination  of  the  diagonal  values  of  matrix  S|_ 
(recall  that  the  diagonal  elements  of  Sl  are  real-valued,  non¬ 
negative,  bounded  by  unity  and  zero,  and  are  arranged  in  order  of 
decreasing  magnitude)  .  As  indicated  earlier,  the  innovations 
model  identified  by  the  algorithm  is  in  balanced  coordinates 

(Moore,  1981) ,  and  the  steady-state  correlation  matrices  of  the 
state  of  both,  the  forward  (II)  and  backward  (nb)  innovations  models 

are  equal  to  matrix  S(_.  In  a  system  representation  in  balanced 
coordinates  the  position  of  a  state  in  the  state  vector  is 
indicative  of  the  importance  of  the  contribution  of  that  state  to 
the  output  correlation  sequence  (the  first  state  is  equal  in 
importance  or  more  important  than  the  second  state;  etc.),  and  the 
magnitude  of  the  corresponding  correlation  matrix  element  is 
representative  of  the  relative  contribution  of  that  state.  Thus, 
an  effective  model  order  selection  approach  is  to  identify  the 
negligible  diagonal  elements  of  matrix  Sj_,  and  select  the  model 

order  to  be  equal  to  the  number  of  non-negligible  diagonal 
elements  of  Sl. 

In  most  situations  involving  a  finite  amount  of  data,  all  the 
diagonal  values  in  matrix  Sl  are  different  from  zero.  This  is  due 

to  the  fact  that  the  subspace  decomposition  is  imperfect  with 
finite  amounts  of  data  because  the  measurement  noise  {w(n)}  corrupts 
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the  past  output  subspace,  and  vice  versa.  In  such  cases,  model 
order  can  be  estimated  by  identifying  jump  discontinuities  in  the 
magnitude  of  the  diagonal  values  of  Sl. 

Other  criteria  can  be  examined  to  estimate  model  order. 
Squaring  the  diagonal  values  of  Sl  emphasizes  discontinuities,  and 

thus  provides  a  good  criterion  also.  The  normalized  running  sum 
of  the  diagonal  values  of  Sl,  and  the  normalized  running  sum  of 
the  squared  diagonal  values  of  Sl,  are  two  additional  criteria  for 
model  order  selection. 


In  the  absence  of  one  or  more  jump  discontinuities,  external 
information  may  be  required,  such  as  prior  knowledge  of  the  system 
being  modeled.  Alternatively,  a  reasonable  model  order  can  be 
selected,  and  various  analyses  can  be  carried  out  to  reduce  the 
order  of  the  model  taking  advantage  of  the  features  of  a  state 
space  realization  in  balanced  coordinates. 


Other  considerations  for  model  order  determination  arise  in 
the  calculation  of  the  QSVD  for  the  matrix  pair  Rg  and  Rq.  It 

turns  out  that  incorrect  determination  of  the  rank  of  a  certain 
matrix  in  the  generation  of  the  QSVD  of  the  matrix  pair  Rg  and  Rq 

can  lead  to  major  difficulties  in  the  determination  of  model 
order.  As  in  Appendix  B,  let  a  2JLxJL  matrix  B  denote  the 
concatenation  of  matrices  Rg  and  Rq  as  (relevant  equations  from 

Appendix  B  are  repeated  here  for  simplicity) 


(3-62) 


Application  of  the  QSVD  to  the  matrix  pair  Rg  and  Rc  leads  to  the 
decomposition,  expressed  for  the  general  case  where  rank(B)<JL, 
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(3-63a) 


B 


UL[SL  0JUL.r(B) 


L  vl[tl 


'JUL-r(B) 


]  J 


vT 


(3-63b) 


Bh  -  Y, 


L'-'JL*  r(B),JLJ 


U 


H 


L  UJL-  r(B),JLJ 


/H 


where  r(B)  denotes  rank(B).  As  indicated  in  Section  3.1,  the  column 
space  of  matrix  Rg  is  equal  to  the  column  space  of  0L,  and  as  is 
well  known,  the  dimension  of  the  column  space  of  0\_  is  equal  to 
the  model  order,  N  .  Therefore,  the  decomposition  of  Rg  in 

Equation  (3-63)  is  indicative  of  model  order.  In  fact,  model 
order  information  is  included  in  matrix  SL,  since  matrix  YL  is  non¬ 
singular  and  matrix  UL  is  unitary.  This  is  another  interpretation 
of  the  order-determining  properties  of  matrix  SL. 

An  important  result  associated  with  the  decomposition  (3-63) 
is  Grassmann's  dimension  theorem.  This  theorem  can  be  stated  as 

(3-64)  dim(range(Rg)  O  range(Rc)]  =  rank(R0)  +  rank(Rc)  -  rank(B) 


Another  relevant  result  follows  directly  from  first  principles, 
(3-65)  rank(B)  £  max[rank(RB),  rank(Rc)] 


Joint  consideration  of  Equations  (3-64)  and  (3-65)  indicates  that 
the  model  order,  N,  satisfies  the  following  constraint, 

(3-66)  N  £  rank(B)  £  JL 

This  equation  implies  that  under-estimation  of  the  rank  of  matrix 
B  in  the  process  of  generating  the  QSVD  for  the  matrix  pair  R0  and 
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Rc  forces  the  dimension  of  the  model  to  be  a  small  value,  whereas 
over-estimation  of  the  rank  of  matrix  B  drives  the  dimension  of 
the  model  towards  large  values.  In  the  approach  adopted  herein 
the  second  condition  is  preferable  because  the  innovations  model 
parameters  are  estimated  in  balanced  coordinates,  and  model  order 
reduction  is  robust  and  effective  in  such  cases,  as  discussed 

above.  Therefore,  in  the  generation  of  the  QSVD  for  the  matrix 
pair  R0  and  Rc  it  is  preferred  to  select  the  rank  of  matrix  B  (or 

the  respective  matrices  in  the  other  QSVDs)  to  be  as  large  as 
possible.  In  most  practical  cases  where  noise  is  present,  the 
rank  is  likely  to  assume  its  maximum  value,  JL. 


4 . 0  INNOVATIONS  SEQUENCE  GENERATION 


In  the  approach  pursued  in  this  program,  an  unknown  system  of 
the  form  (2-2)  is  modeled  as  an  innovations  representation  (2-29) . 
Thus,  once  the  innovations  model  parameters  have  been  identified, 
an  optimal  Kalman  filter  can  be  configured  to  generate  the 
innovations  sequence,  {fi(n)},  for  the  likelihood  ratio  calculations. 

The  approach  described  in  this  section  is  applied  to  the 
observation  data  under  each  of  the  two  hypotheses  . 

Any  one  of  several  equivalent  Kalman  filter  formulations  can 
be  applied  to  generate  the  innovations  sequence.  However,  the 
one-step  predictor  formulation  offers  significant  advantages  in 
the  context  of  the  intended  application  (Anderson  and  Moore, 
1979) .  Specifically,  the  one-step  predictor  formulation  generates 
the  innovations  sequence  and  the  filter  state  update  with  a  simple 
structure  in  ■ the  case  where  the  input  and  output  noises  are 
correlated  (S*[0]  in  Equation  (2-5a)  )  ,  and  thus  imposes  less  real¬ 
time  computational  requirements  than  other  formulations.  Also, 
the  model  identification  algorithm  generates  the  parameters  for 
the  innovations  model.  Thus,  the  one-step  predictor  formulation 
is  adopted  in  this  work.  Strictly  speaking,  the  terminology  "one- 
step  predictor"  should  be  used  hereafter,  but  use  of  the  term 
"Kalman  filter"  is  accepted  universally.  Both  terms  are  used 
herein . 

The  steady-state  one-step  predictor  formulation  for  the 
innovations  model  (2-29)  is  a’  linear,  time-invariant  system 
described  by  the  following  equations: 

(4-ia)  a(n+1|n)  ■  Ffl(n|n-1)  +  K£(n)  n  £  n0 

(4-lb)  £(n) -x(n) -2(n|n-1)  «&(n)  -  HHa(n|n-1)  n  £  nQ 
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( 4-ic)  a(n0|n0*l)«Q 


Here  a(n+1|n)  is  the  estimate  of  the  innovations  model  state  vector 
at  time  n+1  based  on  observation  data  up  to  time  n,  £(n|n-1)  is  the 
estimate  of  the  observation  vector  at  time  n  based  on  observation 
data  up  to  time  n-1,  and  fi(n)  is  the  innovations  associated  with  the 
observation  &(n).  Matrix  K  is  the  steady-state  filter  gain  matrix. 
The  filter  initial  condition  is  set  equal  to  zero  because  the 
innovations  model  initial  condition  is  zero,  Equation  (2-2 9c) .  A 
block  diagram  of  the  Kalman  filter  is  presented  in  Figure  4-1, 
displaying  the  channel  output  vector  as  input,  and  the  innovations 
sequence  vector  as  output . 


Figure  4-1.  Kalman  filter  block  diagram,  emphasizing  the 
innovations  sequence  generation  filter  function. 


The  steady-state  filter  is  an  approximation  to  the  optimal 
time-varying  filter.  If  the  channel  output  process  is  in  steady- 
state,  this  approximation  provides  acceptable  performance. 
Additionally,  the  steady-state  filter  provides  a  significant 
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reduction  in  the  real-time  computational  requirements  over  the 
time-varying  filter.  In  the  cases  where  the  channel  output 
process  is  not  in  steady-state,  filter  performance  is  suboptimal, 
and  the  degree  of  loss  of  optimality  needs  to  be  ascertained. 
Such  a  determination  is  a  subtask  for  Phase  II.  A  related  issue 
involves  filter  initialization  transient  effects.  Since  the 
steady-state  filter  gain  is  used,  it  may  be  necessary  to  neglect 
the  first  Nj  filter  outputs  for  each  data  batch.  Determination  of 
the  value  of  N|  can  be  carried  out  via  analysis  and  simulation,  and 
is  another  subtask  for  Phase  II. 

Anderson  and  Moore  (1979)  show  that  the  filter  estimation 
error  for  an  innovations  model  is  zero  at  all  times.  That  is, 

(4-2)  fl(n+1|n)  »a(n+1) 

Correspondingly,  the  filter  estimation  error  correlation  matrix  is 
zero  also.  This  can  be  inferred  from  the  parallelism  between  the 
innovations  model  (2-29)  and  the  filter  representation  (4-1)  . 
Thus,  knowledge  of  the  filter  implies  knowledge  of  the  innovations 
model,  and  viceversa. 
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5 . 0  LIKELIHOOD  RATIO  DETECTION 


A  detection  methodology  for  complex-valued  multichannel 
Gaussian  processes  has  been  developed  by  Michels  (1991)  in  the 
context  of  innovations-based  detection.  This  approach  has  been 
generalized  recently  to  include  a  class  of  non-Gaussian  processes 
known  as  spherically-invariant  random  processes  (SIRPs)  and  using 
linear  estimators  (Rangaswamy,  Weiner,  and  Michels,  1993)  . 
Michels'  methodology  can  be  applied  directly  to  the  innovations 
sequence  generated  by  the  approach  formulated  herein.  For 
brevity,  only  the  likelihood  ratio  equation  is  presented  here. 

As  discussed  in  Section  4.0,  a  Kalman  filter  (one-step 
predictor)  is  determined  for  each  of  the  two  hypotheses  based  on 
processing  the  multichannel  data.  The  model  order  for  the 
alternative  hypothesis  (H^  filter  is  chosen  to  be  larger  than  the 
model  order  for  the  null  hypothesis  (Hq)  filter.  For  each 
hypothesis  filter,  denote  the  innovations  sequence,  Equation  (4- 
lb) ,  as 

(5-1)  e<n|Hi)  —  x(n)  - 2(n|n-1  ;Ht)  —  x(n)  -  HHa(n|n-1  ;H()  i  - 0, 1 


The  steady-state  correlation  matrix  of  the  innovations  is  denoted 
as  Q(H|),  and  is  defined  in  Equation  (3-60). 


Let  0(H(j,Hi)  denote  the  multichannel  likelihood  ratio  as 
defined  by  Michels  (1991)  for  the  Gaussian  signal  case.  Then,  the 
log-likelihood  ratio  (LLR)  can  be  expressed  as  follows, 


Nr 


:5-2>  ln[e(H„,H,)]  .  £ 

n«n. 


Ini 


I  QjWfl)  I 


n(H,)|J 


+  £H(n|H0)  n  ^Ho)  £<n|H0) 
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The  LLR  is  compared  to  a  threshold,  T,  which  is  calculated 
adaptively  to  maintain  a  constant  false  alarm  rate  (CFAR) , 


(5-3) 


lr[e(H0,H,)] 


£  T 
<  T 


select  H1 
select  H0 


A  candidate  CFAR  approach  with  demonstrated  good  performance 
calculates  the  median  of  a  set  of  the  LLR  values  from  a  number  of 
adjacent  range  cells  (at  the  same  azimuth)  on  both  sides  of  the 
cell  in  question,  and  scales  the  calculated  median  value  by  a 
pre-determined  constant  to  provide  the  desired  false  alarm  rate 
(Metford  and  Haykin,  1985)  . 

The  LLR  expression  has  to  be  modified  if  optimal  time-varying 
filters  are  used  instead  of  the  steady-state  filters.  In  such 
cases  the  modification  is  straightforward,  and  involves  replacing 
the  steady-state  correlation  matrices  of  the  two  innovations  by 
their  time-varying  values. 

Alternative  expressions  for  the  log-likelihood  ratio  can  be 
generated  based  on  factorization  of  the  innovations  correlation 
matrix  and  spatial  whitening  of  the  innovations  process.  This 
includes  Cholesky  factorization,  LDU  decomposition,  and  SVD .  The 
first  two  techniques  have  been  described  by  Michels  (1991),  and 
lead  to  simplified  LLR  expressions.  The  SVD  technique  is  derived 
here . 


Consider  the  steady-state  innovations  correlation  matrices 
for  each  of  the  two  hypotheses  and  carry  out  an  SVD  on  each 
correlation  matrix.  This  results  in  the  following  decompositions: 

(5-4)  a(Hj)  .  VjljVf  i-0,  1 
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where  matrix  V|  is  a  JxJ  unitary  matrix,  and  £j  is  a  diagonal  matrix 
with  real-valued,  positive  elements  arranged  along  the  diagonal  in 
decreasing  order  of  magnitude  (it  is  assumed  herein  that  the 
correlation  matrix  of  the  innovations  sequence  has  full  rank) . 


That  is, 


L  0  0  ...  cj 


(5-5b)  oj*  i  ojj  *  . . .  *  Oy  >  0 

Since  matrix  V|  is  unitary,  the 
of  Q(H|)  are  obtained  easily  as 

(5-6)  OT1  (H|)  -  V, Ij1 

(5-7)  |Q(H|)|  -  fl  oj 

k-1 


i  ■  0, 1 

1-0.1 

determinant  and  inverse  functions 

i-0. 1 

i  m  0,  1 


Now  make  a  linear  transformation  on  the  innovations  sequence  using 
the  unitary  matrix  Vj,  to  obtain 

(5-8)  Jt(n|H,)  -  vfc(n|H|)  1-0,1 

The  transformed  innovations  sequence,  ¥.(n|Hj),  is  uncorrelated 
spatially  as  well  as  temporally  (recall  that  fc(n|Hj)  is  uncorrelated 
temporally),  with  correlation  matrix  £j.  Transformation  of  a  J- 
dimensional  vector  by  a  unitary  matrix  rotates  the  vector  in  the  J- 
dimensional  space,  but  does  not  alter  its  magnitude.  Thus,  the 
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spatial  whitening  transformation  does  not  alter  the  variance  of 
the  elements  of  the  innovations  vector. 


Substituting  Equations  (5-4)  through  (5-8)  into  Equation  (5- 
2)  results  in  the  following  LLR  expression 


(5-9) 


lr[e(H0.M1)]  -St 

n-%  k-1 


In 

.  1  vk(n|H0 )  |2  |vk(n|H,)|2' 

2  2 
°0k  °1k  J 

where  Vk(n|H|)  denotes  the  Kth  element  of  y(n|Hj).  This  LLR  is  of  the 

same  form  as  the  LLR  derived  by  Michels  (1991)  for  spatial 
whitening  of  the  innovations  using  an  LDU  decomposition. 
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6 . 0  SOFTWARl  SIMULATION 


The  identification  and  filtering  algorithms  described  in  the 
preceding  sections  have  been  programmed  in  FORTRAN  77  for  Apple 
Macintosh  processors.  Support  software  for  the  validation  and 
execution  of  the  routines  has  been  generated  also.  The  support 
software  includes  signal  generation  routines/  auxiliary  routines 
for  validation^  and  code  for  miscellaneous  calculations.  The 
identification  algorithm  makes  use  of  the  singular  value  and  QR 
decompositions.  Subroutines  that  implement  these  matrix 
operations  for  complex-valued  matrices  were  obtained  from  versions 
of  the  LINPACK  software  package  (Dongarra  et  al.,  1979) .  Separate 
code  was  written  and  exercised  to  validate  the  LINPACK  routines 
before  incorporation  into  the  main  algorithm  code.  The  signal 
generation  code  uses  a  Gaussian  random  number  generator  obtained 
from  the  text  by  Press  et  al.  (1989)  .  Sample  realizations 
generated  by  this  code  were  tested  for  whiteness  and  gaussianity. 

*  .  1  flaftirtri _ Validation 

Code  validation  was  carried  out  in  two  steps .  First,  all 
subroutines  and  select  segments  of  code  were  validated 
individually.  Second,  the  complete  package  was  validated  using 
examples  generated  for  that  purpose.  The  examples  consisted  of 
system  models  with  a  simple  structure  so  that  the  computer  output 
could  be  predicted.  Both  real-valued  and  complex-valued  examples 
were  generated.  One  particular  example  used  is  the  second-order 
system  defined  by  the  following  matrix  parameters  (for  a  system 
model  of  the  form  (2-2)): 


L  ^21  ^22  -I 
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Hh  -  G  -  Dh  -  Q  -  C  -  l2 

This  model  was  used  to  generate  a  random  vector  sequence  to 
validate  various  aspects  of  the  software.  For  example,  defining 
matrix  F  with  ■  f12  ■  I22  m  0  and  ^21*^  generates  an  output  vector 

sequence  that  consists  of  white  noise  in  each  output  channel,  but 
the  two  channels  are  correlated  from  one  instant  to  the  next  (the 
correlation  is  due  to  the  coupling  induced  by  the  non-zero  (2,1) 
element  of  F)  .  The  output  of  the  identification  program  should 

indicate  a  first-order  model,  with  the  first  diagonal  element  of 
matrix  Sl  approximately  equal  to  0.7071,  and  low  values  for  the 

remaining  diagonal  elements.  This  was  the  result  obtained. 
Complex-valued  test  cases  using  this  sample  model  were  generated 
by  letting  F  be  a  diagonal  matrix  with  the  desired  complex-valued 
poles  along  the  diagonal. 

During  validation  and  testing  it  was  discovered  that  system 
poles  along  the  real  axis  are  more  difficult  to  estimate,  and  that 
Equation  (3-49)  can  produce  biased  results  (over-estimation  of  the 
system  poles)  in  some  cases.  This  has  been  observed  also  by  the 
research  team  at  the  Catholic  University  of  Leuven,  and  one 
approach  to  mitigate  this  consists  of  utilizing  alternative,  more 
complex  algorithms  for  the  estimation  of  matrix  F.  Identification 
and  detailed  evaluation  of  these  alternatives  is  an  activity  for 
Phase  II  of  this  program. 

6 . 2  ftaalyafla _ and  Simulation  Results 

The  software  has  been  exercised  with  cases  generated  using 
multichannel  AR  models  provided  by  the  program  monitor  at  RL,  Dr. 
James  H.  Michels.  These  cases  consist  of  signal  only,  clutter 
only,  signal  +  noise,  clutter  +  noise,  and  signal  +  clutter  + 
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noise.  In  all  casas  tha  signal,  clutter,  and  noise  processes  are 
statistically  independent  of  each  other. 

Signal  Aft  Mftdol 


The  signal  model  is  a  complex-valued,  two-input,  two-output 
AR  model  of  order  2  with  the  following  matrix  parameters, 

&(n)  -  *  Aj(1)&(n-1)  -  Ay(2)fc(n*2)  +  u,(n) 


A?(1) 


1.6290  *j  1.4241  x10'7 
1 . 3733x1 O'®  +  j  3.8202x1 O'13 


1.3733x10'®  +  j  3.8202x1  O'13 
1.6290  -  j  1.4241x1 0*7 


0.80996  -  j  1.41 62x1  O'7  1.5259x10'®  -  j  9.0949x1  O'13' 

1.5259x10'® -j  9.0949x1  O'13  0.80996  -  j  1 .4162x10‘7 


The  input  to  the  signal  AR  recursion,  (U$(n)},  is  a  zero-mean,  unit 
variance  white  noise  sequence  with  a  spatial  correlation  structure 
defined  as 


Q. 


0.13038  0.12907 

0.12907  0.13038 


This  two-input,  two-output  AR  model  corresponds  to  a  fourth-order 
state  space  model  in  an  innovations  representation  (as  described 
in  Appendix  A)  ,  with  poles  at  the  following  locations  in  the 
complex  z-plane: 


True  Signal  Model  Poles:  -0.81451  ±  j  0.38282 

-0.81449  ±  j  0.38281 


This  AR  system  was  defined  by  Michels  to  have  a  very  high  channel- 
to-channel  correlation  (**0.99),  which  indicates  that  a  lower-order 
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modal  could  represent  the  signal  information.  Specifically,  a 
second-order  state  space  model  can  represent  the  signal 
information  well.  Notice  that  the  pole  locations  are  almost 
repeated  roots,  which  indicates  that  the  two  channels  are  almost 
identical.  Therefore,  given  a  high-level  of  channel-to-channel 
correlation,  a  reduced-order  model  should  perform  adequately. 

The  AR  process  {j£*(n)}  is  corrupted  by  a  zero-mean,  unit- 
variance  white  noise  sequence  {lflt(n)}  to  give  the  noise-corrupted 
channel  output  sequence  as 

n)  •  JU(n)  +  jfld(n) 

For  this  noise  model  and  the  signal  model  given  above,  the  signal- 
to-noise  ratio  (SNR)  is  approximately  3  dB. 

Consider  the  problem  of  representing  the  AR  signal  in 
additive  white  noise  with  a  state  space  model.  The  channel  output 
noise,  {)flt(n)},  alters  the  parameters  of  the  state  space  model 
designed  for  the  AR  signal  {^(n)}  only,  but  {**(")}  can  be  represented 
as  the  output  of  a  state  space  model.  That  is,  (Z$(n)}  is 

represented  as  the  output  of  an  innovations  model,  but  the  model 
for  {*e(n)>,  which  includes  the  additive  noise  {W.(h)},  is  not  an 
innovations  model  (there  is  an  innovations  model  for  {&g(n)},  but  it 
is  different  from  the  innovations  model  for  tes(n)}).  This  is  a 
manifestation  of  the  well-known  fact  that  an  AR  process  corrupted 
by  additive  output  white  noise  is  no  longer  an  AR  process.  In 
contrast,  the  state  space  model  remains  a  valid  representation  of 
the  signal  even  after  the  addition  of  a  new  noise  source.  This  is 
one  of  the  reasons  why  algorithms  developed  based  on  state  space 
models  are  more  robust  than  algorithms  based  on  time  series 
models . 
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Clutter  AR  Model 


The  clutter  model  is  a  complex-valued,  two-input,  two-output 
AR  model  of  order  2  with  the  following  matrix  parameters, 

¥c(n)  -  -  Ac(1)jk(n-1)  -  Ac(2)&;(n-2)  +  y^n) 


Ac(1)  - 


-1 .0430 

0.0 


0.0 

-1.0430 


*8(2) 


0.4900  0.0 

.  0.0  0.4900. 


The  input  to  the  clutter  AR  recursion,  {lic(n)}f  is  a  zero-mean,  unit 
variance  white  noise  sequence  with  a  spatial  correlation  structure 
defined  as 


Qc 


'1.5502 

0.0 


0.0 

1.5502. 


This  two-input,  two-output  AR  model  corresponds  to  a  fourth-order 
state  space  model  in  an  innovations  representation  (see  Appendix 
A),  with  poles  at  the  following  locations  in  the  complex  Z-plane : 

True  Clutter  Model  Poles :  0.5215  ±  j  0.4669 

0.5215  ±  j  0.4669 

The  clutter  AR  coefficient  values,  the  noise  covariance  values, 
and  the  diagonal  structure  of  this  AR  system  indicate  that  the  two 
channels  are  uncorrelated.  Thus,  a  fourth-order  state  space  model 
can  represent  the  clutter  information  well.  Notice  that  the  pole 
locations  are  repeated  roots,  which  indicates  that  the  two 
channels  are  identical. 
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The  clutter  AR  process  {Yc<n)}  is  corrupted  by  the  zero-mean, 
unit-variance  white  noise  sequence  {W.(n)}  to  give  the  noise- 
corrupted  channel  output  sequence  as 

Xc(n)  -  yc(n)  +  Jfld(n) 


For  this  noise  model  and  clutter  model  the  clutter-to-noise  ratio 
(CNR)  is  approximately  6  dB. 

Selected  Simulation  Results 

The  identification  and  filtering  software  was  exercised  with 
the  sequence  Us(n)}  •  Calculated  values  of  the  model  order  criterion 
parameters,  the  diagonal  elements  of  matrix  S|_,  indicate  that  a 
second-order  state  space  model  is  a  good  approximation  to  this 
system,  as  expected.  Specifically,  the  values  of  the  diagonal 
elements  of  matrix  Sl  are:  {0.9804,  0.9776,  0.6712,  0.6360, 

0.1215,  0.0912},  and  the  set  of  the  square  of  these  values  is: 
{0.9612,  0.9557,  0.4505,  0.4045,  0.0148,  0.0083}.  Examination  of 
these  two  sets  indicates  that  a  second-order  model  is  a  good  fit 
to  the  data.  The  next  reasonable  model  order  selection  is  four. 

Based  on  the  above  discussions,  model  order  2  was  selected 
for  the  AR  signal  in  white  noise.  Figure  6-1  is  a  plot  of  the 
real  part  of  the  first  element  of  a  single  realization  of  the 
innovations  vector  process,  (£i(n)},  generated  using  a  filter  of 

order  2  (all  plots  herein  are  for  single-realization  cases) .  The 
filter  parameters  were  identified  using  2,214  output  sequence 
vectors  (corresponding  to  L  =  8  and  M=  2,200  in  Equations  (3-21) 
and  (3-22) ) .  Only  500  points  are  shown  in  the  figure,  but  these 
points  are  representative  of  the  sample  process.  The  innovations 
sequence  appears  to  be  unbiased,  with  a  calculated  sample  mean  of 
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e(n) 


0.0103  +  j  0.0074' 
.  0.0308  -  j  0.0050 . 


Notice  also  the  high  degree  of  "whiteness"  exhibited  by  the 
innovations.  The  imaginary  part  of  the  innovations  behaves 
similarly.  The  autocorrelation  function  of  the  sequence  in  Figure 
6-1  was  estimated,  and  is  shown  in  Figure  6-2,  This  figure 
clearly  indicates  the  random  (white)  nature  of  the  innovations,  as 
expected.  The  zero-lag  innovations  correlation  matrix  identified 
by  the  software  using  Equation  (3-60)  is 

0  [  1.5337  0.5908  +  j  0.0453 
■  *  [  0.5908  -  j  0.0453  1 .5832  J 

This  agrees  very  well  with  the  sample  correlation  value  of  1.536 
indicated  in  Figure  6.2.  Several  simulation  runs  were  made  using 
multiple  sample  realizations  of  the  same  length  and  filter  order 
two.  A  plot,  of  the  innovations  correlation  averaged  over  ten 
realizations  looks  very  similar  to  Figure  6-2. 


Figure  6-1.  Real  part  of  the  first  element  of  innovations  vector 

for  the  case  of  signal  plus  noise. 
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first  element  of  innovations  covariance  for  case  of  signal  plus  noise 


10 

time  index,  n 


Figure  6-2.  Real  part  of  the  auto-correlation  function  of  the 
first  element  of  the  innovations  sequence  vector  for  the  case  of 

signal  plus  noise. 


Identification  algorithm  performance  can  be  assessed  by 
examining  the  roots  of  the  identified  innovations  model  system 
matrix,  F.  The  scatter  plots  in  Figure  6-3,  which  correspond  to 
results  obtained  for  ten  distinct  realizations,  illustrate  the 
parameter  identification  capability  of  the  algorithm.  These 
scatter  plots  show  the  ten  identified  root  pairs,  all  in  close 
proximity  to  the  true  roots  given  above.  The  dashed  lines  in  each 
plot  intersect  at  the  center  of  each  plot  (-0.81  -  j  0.38  and  - 
0.81  +  j  0.38,  respectively),  and  the  centers  are  very  close  to 
the  true  root  values  (-0.8145  ±  j  0.3828).  All  the  identified 
roots  are  at  a  distance  less  than  3%  of  the  true  values,  and  most 
are  much  closer  than  that . 
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Figure  6-3.  Scatter  plot  of  real  and  imaginary  parts  of 
identified  model  poles  for  ten  distinct  realizations  of  signal 

plus  noise. 


The  software  was  used  also  to  model  and  analyze  the  clutter 
plus  noise  sequence,  tSc(n)}.  For  this  case  at  a  CNR  of  20  dB,  the 

model  order  criterion  parameters,  the  diagonal  elements  of  the 
matrix  SL  are:  {0.8187,  0.7919,  0.3027,  0.2748,  0.1340,  0.1112); 

and  the  set  of  the  square  of  these  values  is:  {0.6703,  0.6271, 

0.0916,  0.0755,  0.0180,  0.0124}.  This  information,  together  with 
knowledge  of  the  lack  of  channel  correlation,  indicates  that  a 
fourth-order  state-space  model  is  a  good  approximation  to  this 
system.  Without  prior  knowledge  regarding  channel  correlation, 
measures  such  as  the  percentage  incremental  power  attributable  to 
each  additional  singular  value  may  result  in  improved  model  order 
estimates . 


Model  order  four  was  selected  for  state  space  representation 
of  the  clutter  AR  process  in  additive  white  noise.  Plots  of  the 
real  and  imaginary  parts  of  the  first  element  of  a  single 
realization  of  the  innovations  vector  process,  {£i(n)},  are  presented 
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in  Figure  6.4.  These  results  were  generated  using  a  fourth-order 
filter  for  a  case  with  6  dB  CNR.  The  filter  parameters  were 
identified  using  2,214  output  sequence  vectors  (corresponding  to  L 
»  8  and  2,200  in  Equations  (3-21)  and  (3-22)),  as  in  the  signal 
+  noise  case.  Only  500  points  are  shown  Figure  6.4,  but  these 
points  are  representative  of  the  sample  process.  Both  components 
(real  and  imaginary)  of  the  sequence  {e^  (n)}  are  unbiased,  as 

indicated  by  the  sample  mean, 


I(n) 


-0.0101  -j  0.0429 
0.0538  +  j  0.0076 


An  estimate  of  the  real  and  imaginary  parts  of  the  sample 
autocorrelation  function  of  {^(n)}  of  Figure  6-4  is  given  in  Figure 
6-5.  The  real  part  has  an  impulse  at  lag  n = 0  and  is  close  to  zero 
everywhere  else,  which  is  representative  of  a  white  innovations. 
The  imaginary  part  exhibits  low-amplitude  oscillations  about  zero, 
as  expected  of  a  white  innovations .  The  zero-lag  innovations 
correlation  matrix  estimated  using  Equation  (3-60)  is 


a 


3.1114  0.0594  +j  0.0248' 

0.0594  -  j  0.0248  3.3048 


Element  (1,1)  agrees  with  the  sample  correlation  value  of  3.106  + 
j  0.0  indicated  in  Figure  6-5.  The  behaviour  of  {£2(0)}  is  similar. 

Figure  6-6  presents  scatter  plots  of  the  poles  of  the  fourth- 
order  system  for  ten  realizations.  The  roots  are  clustered  about 
the  values  of  the  true  repeated  roots  (0.5215  ±  j  0.4669),  which 
are  indicated  by  the  intersections  of  the  dashed  lines .  The 
largest  root  estimation  error  is  less  than  8.2%.  This  error  is 
larger  than  the  worst  error  in  the  signal  plus  noise  case,  and  is 
due  to  the  greater  difficulty  in  estimating  faster  modes. 
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imaginary  component  of  innovations 


first  element  of  innovations  vector  for  case  of  clutter  plus  noise 


first  element  of  innovations  vector  for  case  of  clutter  plus  noise 


Figure  6-4.  Real  and  imaginary  parts  of  the  first  element  of  the 
innovations  sequence  vector  for  the  case  of  clutter  plus  noise 

(CNR  -  6  dB  conditions) . 
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imaginary  component  of  covariance  real  component  of  covariance 


first  element  of  innovations  covariance  for  null  hypothesis  data 
using  null  hypothesis  filter  (order  4,  CNR=6dB) _ 
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Figure  6-5.  Real  and  imaginary  parts  of  the  auto-correlation 
function  of  the  first  element  of  the  innovations  sequence  vector 
for  the  case  of  clutter  plus  noise  (CNR  -  6  dB) . 
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Figure  6-6.  Scatter  plot  of  real  and  imaginary  parts  of 
identified  model  poles  for  ten  distinct  realizations  of  clutter 

plus  noise. 
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Similar  biased  behaviour  has  been  observed  in  other 
estimation  results  (Michels,  1992b),  as  well  as  in  detection 
performance  results  (Michels,  1992a)  obtained  using  time  series 
(AR)  models.  For  the  state-space  approach  pursued  herein, 
unbiased  estimates  with  reduced  variance  can  be  obtained  by 
averaging  several  individual  estimates  and/or  by  increasing  the 
complexity  of  the  estimation  algorithm. 

Various  simulations  were  carried  out  to  obtain  a  first-order 
assessment  of  the  discrimination  capability  of  the  innovations- 
based  methodology  using  the  SSC  algorithm.  One  set  of  simulations 
involved  designing  a  Kalman  filter  for  each  hypothesis,  processing 
data  corresponding  to  each  of  the  two  hypotheses  using  both 
filters,  and  analyzing  the  resulting  four  filter  output  sequences 
(two  filters,  and  each  filter  processes  data  sets  corresponding  to 
each  of  the  two  hypotheses).  These  results  are  presented  next. 
As  before,  all  plots  correspond  to  single-realization  cases. 

Consider  first  the  case  of  processing  data  from  each  of  the 
two  hypotheses  using  a  null  hypothesis  filter,  corresponding  to 
clutter  +  noise  only.  For  this  case  the  filter  order  is  four,  as 
mentioned  earlier  in  the  clutter  plus  noise  model  discussion. 
Results  are  presented  herein  for  two  sets  of  conditions:  (a)  SNR  = 
3  dB  and  CNR  -  6  dB;  and  (b)  SNR  -  3  dB  and  CNR  -  20  dB .  For  each 
set  of  conditions  the  procedure  described  next  was  followed. 

•  A  realization  of  the  clutter  +  noise  process  of  duration 

M  ■  2,200  was  generated  and  processed  to  design  a 

fourth-order  Kalman  filter.  The  resulting  filter  is  the 
filter  for  the  null  hypothesis  (signal  not  present) . 

•  The  null  hypothesis  filter  was  applied  to  a  clutter  + 

noise  process  sequence  of  duration  2,200,  and  the 
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sample  correlation  matrix  sequence  of  the  filter  output 
sequence  was  calculated.  The  real  and  imaginary  parts 
of  the  (1,1)  element  of  the  resulting  sample  correlation 
matrix  sequence  are  plotted  in  Figure  6-5  for  CNR  -  6  dB 
conditions,  and  Figure  6-7  for  CNR  -  20  dB  conditions. 

Both  sets  of  figures  are  representative  of  the  auto¬ 
correlation  of  a  white  innovations  sequence,  as  expected 
(both  sets  of  figures  show  low-level  energy  content  at 
the  higher  lags) . 

•  The  null  hypothesis  filter  was  applied  to  a  combined 
signal  +  clutter  +  noise  process  sequence  (alternative 
hypothesis  case)  of  duration  2,200,  and  the  sample 

correlation  matrix  sequence  of  the  filter  output 
sequence  was  calculated.  In  this  case,  however,  the 
sequence  is  not  a  true  innovations  sequence  because  the 
filter  is  not  optimal  for  this  process.  The  real  and 
imaginary  parts  of  the  (1,1)  element  of  the  resulting 
sample  correlation  matrix  sequence  are  plotted  in  Figure 
6-8  for  CNR  -  6  dB  conditions,  and  Figure  6-9  for  CNR  - 
20  dB  conditions.  Both  of  these  figures  show  a  marked 
deviation  from  the  expected  auto-correlation  for  a  white 
innovations  sequence. 

In  the  discussions  and  results  presented  above  the  (2,2)  element 
of  the  sample  correlation  matrix  is  not  referred  to.  This  is  due 
to  the  fact  that  its  behaviour  is  very  similar  to  the  behaviour  of 
the  (1,1)  element. 

In  continuation  of  the  first-order  assessment  of  the 
discrimination  capability  of  the  SSC  approach,  consider  now  the 
case  of  processing  data  from  each  of  the  two  hypotheses  using  an 
alternative  hypothesis  filter,  corresponding  to  the  combined 


process  of  signal  +  clutter  +  noise.  Since  the  signal  and  clutter 
are  uncorrelated  m  this  set  of  examples,  a  sixth-order  state 
space  model  is  required  tor  the  combined  process.  As  before, 
results  are  presented  for  two  sets  of  conditions:  (a)  SNR  ■  3  dB 
and  CNR  ■  6  dB;  and  (b)  SNR  -  3  dB  and  CNR  «  20  dB.  For  each  set 
of  conditions  the  procedure  described  next  was  followed  (all  plots 
are  for  single-realization  cases) . 

•  A  realization  of  the  combined  signal  +  clutter  +  noise 

process  of  duration  M  ■  2,200  was  generated  and 
processed  to  design  a  sixth-order  Kalman  filter.  The 

resulting  filter  is  the  filter  for  the  alternative 
hypothesis  (signal  present)  . 

•  The  alternative  hypothesis  filter  was  applied  to  a 

coi.jained  process  sequence  of  duration  M»  2,200,  and  the 

sample  correlation  matrix  sequence  of  the  filter  output 
sequence  was  calculated.  The  real  part  of  the  (1,1) 
element  of  the  resulting  sample  correlation  matrix 
sequence  is  plotted  in  Figure  6-10  for  CNR  -  6  dB 

conditions,  and  Figure  6-11  for  CNR  =*  20  dB  conditions. 

As  attested  in  both  figures,  the  correlation  sequences 
correspond  to  white  innovations  sequences,  as  expected 
(both  figures  show  low-level  energy  content  at  the 
higher  lags) . 

•  The  alternative  hypothesis  filter  was  applied  to  a 

clutter  +  noise  process  sequence  (null  hypothesis  case) 
of  duration  M*  2,200,  and  the  sample  correlation  matrix 
sequence  of  the  filter  output  was  calculated.  In  this 
case,  however,  the  sequence  is  not  a  true  innovations 

sequence  because  the  filter  is  not  optimal  for  this 
process.  The  real  part  of  the  (1,1)  element  of  the 


resulting  sample  correlation  matrix  sequence  is  plotted 
in  Figure  6-12  for  CNR  ■  6  dB  conditions#  and  Figure  6- 
13  for  CNR  »  20  dB  conditions.  The  correlation  sequence 
in  each  of  the  figures  corresponds  to  a  colored  process# 
and  not  to  a  white  innovations  sequence.  This  is  the 
expected  result . 

Figures  6-10  through  6-33  do  not  include  the  imaginary  part  of  the 
sample  correlation  sequence  because  it  is  similar  to  the  imaginary 
part  of  the  sample  correlation  sequence  presented  in  the  preceding 
figures.  For  this  case  also  the  behaviour  of  the  (2# 2)  element  is 
very  similar  to  that  of  the  (1#1)  element  of  the  sample 
correlation  matrix  sequence 

These  results  indicate  that  the  innovations-based  detection 
methodology  using  the  identification  algorithm  adopted  in  this 
program  can  discriminate  between  data  corresponding  to  each  of  the 
two  hypotheses.  That  is,  a  filter  designed  for  the  alternative 
hypothesis  (signal  +  clutter  +  noise)  generates  a  true  innova*  ons 
sequence  given  a  signal  +  clutter  +  noise  channel  process,  und 
generates  a  colored  output  given  a  clutter  +  noise  channel 
process.  Analogously,  a  filter  designed  for  the  null  hypothesis 
(clutter  +  noise)  generates  a  true  innovations  sequence  given  a 
clutter  +  noise  channel  process,  and  generates  a  colored  output 
given  a  signal  +  clutter  +  noise  channel  process. 
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Figure  6-7 .  Real  and  imaginary  parts  of  the  auto-correlation 
function  of  the  (1,1)  element  of  the  innovations  sequence  vector 
for  the  case  of  null  hypothesis  data  using  the  null  hypothesis 
filter  (CNR  =  20  dB  conditions) . 
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Figure  6-8.  Real  and  imaginary  parts  of  the  auto-correlation 
function  of  the  (1,1)  element  of  the  filter  output  vector  for  the 
case  of  alternative  hypothesis  data  using  the  null  hypothesis 
filter  (CNR  =  6  dB  conditions) . 
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Figure  6-9.  Real  and  imaginary  parts  of  the  auto-correlation 
function  of  the  (1,1)  element  of  the  filter  output  vector  for  the 
case  of  alternative  hypothesis  data  using  the  null  hypothesis 
filter  (CNR  =  20  dB  conditions) . 
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first  element  of  innovations  covariance  for  alternative 
hypothesis  data  using  alternative  hypothesis  filter  (order  6 


time  index,  n 

Figure  6-10.  Real  part  of  the  auto-correlation  function  of  the 
(1,1)  element  of  the  innovations  sequence  vector  for  the  case  of 
alternative  hypothesis  data  using  alternative  hypothesis  filter 

(CNR  -  6  dB  conditions) . 
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Figure  6-11.  Real  part  of  the  auto-correlation  function  of  the 
(1,1)  element  of  the  innovations  sequence  vector  for  the  case  of 
alternative  hypothesis  data  using  alternative  hypothesis  filter 

(CNR  -  20  dB  conditions) . 


80 


c 

Q> 

C 

8L 

E 

8 

73 

2 


first  element  of  innovations  covariance  for  null  hypothesis 
data  using  alternative  hypothesis  filter  (order  6) _ 
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Figure  6-12.  Real  part  of  the  auto-correlation  function  of  the 
(1,1)  element  of  the  filter  output  vector  for  the  case  of  null 
hypothesis  data  using  the  alternative  hypothesis  filter 

(CNR  =  6  dB  conditions) . 
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Figure  6-13.  Real  part  of  the  auto-correlation  function  of  the 
(1,1)  element  of  the  filter  output  vector  for  the  case  of  null 
hypothesis  data  using  the  alternative  hypothesis  filter 

(CNR  =  20  dB  conditions) . 


7 . 0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  work  carried  out  in  this  program  emphasized  the 
development  and  analysis  of  a  state  space  methodology  and 
algorithm  for  the  model-based  multichannel  detection  problem  in 
the  context  of  radar  system  applications.  Application  of  state 
space  techniques  for  multichannel  detection  in  radar  systems  is 
one  novel  aspect  of  the  work  reported  here.  The  state  space  model 
class  is  richer  than  the  time  series  model  class  that  is  used 
often  in  radar  system  applications.  And,  as  demonstrated  in  this 
work,  the  state  space  model  class  can  be  used  to  represent 
effectively  multichannel  radar  signals. 

Another  novel  aspect  of  the  work  is  t-he  utilization  in  the 
detection  methodology  of  a  new  algorithm  developed  by  Van 
Overschee  and  De  Moor  (1993) .  This  algorithm  was  adopted  in  the 
program  for  .multichannel  radar  output  modeling  and  parameter 
identification.  In  the  process,  the  algorithm  was  extended  to  the 
case  of  complex-valued  radar  system  data,  and  an  alternative 
derivation  of  the  algorithm  was  developed  which  is  simpler  and 
easier  to  follow  than  the  one  presented  in  the  forthcoming  paper 
(Van  Overschee  and  De  Moor,  1993) .  The  selected  approach  uses 
channel  output  data  directly  (as  opposed  to  output  correlation 
matrices)  to  estimate  model  parameters.  This  eliminates  the  large 
computational  burden  associated  with  the  generation  of  the  output 
correlation  matrix  sequence,  and  leads  to  reduced  numerical 
precision  (dynamic  range)  requirements.  Furthermore,  in  a 
practical  environment  it  may  be  possible  to  start  processing  the 
data  as  it  is  received.  In  contrast,  techniques  which  require  the 
computation  of  channel  output  correlation  matrices  have  a  built-in 
delay  because  the  calculation  of  every  lag  requires  availability 
of  all  the  channel  output  sequence. 
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The  Van  Overschee-De  Moor  algorithm  belongs  to  a  class  of 
techniques  referred  to  as  subspace  methods.  Subspace  methods  are 
based  on  decomposing  the  vector  space  spanned  by  the  channel 
outputs  into  signal  and  noise  subspaces.  This  decomposition  is 
carried  out  with  robust  numerical  techniques  such  as  the  SVD  and 
the  QR  decomposition.  Thus,  the  algorithm  offers  numerical  and 
performance  advantages  over  other  techniques. 

A  computer  simulation  was  developed  to  validate  the  algorithm 
and  methodology,  and  to  serve  as  a  testbed  for  evaluation  of  the 
algorithm  in  radar  system  applications.  The  simulation  can  be 
exercised  with  internally-generated  sample  multichannel  output 
data,  or  with  externally-provided  data.  Extensive  tests  were 
carried  out  to  validate  the  code. 

Simulation-based  analyses  carried  out  to  date  demonstrate  the 
feasibility  of  the  SSC  state  space  approach  for  multichannel 
identification  and  detection  in  radar  system  applications.  The 
algorithm  has  demonstrated  the  capability  to  discriminate  between 
signal  plus  clutter  plus  noise  and  clutter  plus  noise  in  an 
innovat ions-based  detection  algorithm  formulation  for  the 
multichannel  detection  problem.  Several  cases  have  been  analyzed 
at  various  SNR  and  CNR  levels,  and  in  all  cases  simulated  thus  far 
discrimination  has  been  demonstrated. 

In  the  process  of  completing  the  work  reported  here  several 
areas  have  been  identified  for  further  research  and  development  in 
a  Phase  II  of  this  program.  These  areas  are  summarized  below. 

Processor.  Requirements  Definition 

Determination  of  the  true  potential  of  the  SSC  approach  for 
radar  system  applications  requires  the  establishment  of  a  detailed 
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set  of  requirements  for  various  radar  problems  such  as  space/time 
processing  in  a  radar  array  system,  the  fusion  of  data  from 
multiple  distinct  rk iar  systems,  and  the  fusion  of  multiple  radar 
measurements  (dual-polarization  measurements,  for  example) .  Also, 
detailed  sets  of  requirements  need  to  be  defined  for  other 
application  areas  such  as  hydrological  systems,  seismic  event 
detection,  and  medical  technology. 

Additional  Analyses  and  Detailed  Algorithm  Formulation 

The  analyses  listed  below  are  required  to  generate  a  detailed 
algorithm  definition  for  the  requirements,  and  to  provide  a 
precise  assessment  of  the  SSC  approach  in  the  context  of  the 
requirements . 

•  The  innovations  model  matrix  parameters  F,  T,  and  H  can 
be  estimated  using  different  equations.  Some  of  the 
alternative  equations  can  exhibit  bias  errors,  but  may 
be  simpler  to  implement.  These  alternatives  need  to  be 
evaluated  and  traded. 

•  Model  order  selection  criteria  for  on-line  and  off-line 
decisions  need  to  be  evaluated  and  traded.  This 
includes  the  diagonal  values  { Sj}  of  matrix  Sj_,  the 

square  of  the  {Sj},  and  their  normalized  running  sums. 

•  The  steady-state  Kalman  filter  was  used  in  this  Phase  I 
to  generate  the  innovations  sequence.  Alternatively, 
the  time-varying  Kalman  filter  can  be  used.  The  loss 
in  performance,  if  any,  incurred  by  using  the  steady- 
state  approximation  needs  to  be  evaluated.  A  related 
issue  is  the  duration  of  the  transient  effect  in  the 
case  of  the  steady-state  filter. 
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•  Key  implementation  parameters  need  to  be  established. 

This  includes  the  minimum  required  channel  output 
sequence  duration,  and  the  row  dimension  of  the  Hankel 
output  data  matrix. 

•  Identification  and  detection  performance  should  be 
compared  with  that  of  other  methods.  This  includes 
state  space  methods  that  operate  on  output  correlation 
matrix  data,  and  methods  based  on  AR  models. 

•  A  detailed  algorithmic  approach  to  the  implementation 
of  the  QR  decomposition  and  the  QSVD  needs  to  be 
defined.  in  this  context,  a  new  decomposition  for 
rectangular  matrices  introduced  by  Stewart  (1992) 
should  be  reviewed  for  possible  utilization  in  the  SSC 
approach.  This  decomposition  is  related  to  the  QR 
decomposition  and  to  the  SVD,  and  can  be  updated 
recursively  in  a  simple  manner.  The  recursive  feature 
is  attractive  for  reducing  the  computational  load. 

Once  these  technical  issues  are  addressed,  a  detailed  architecture 
design  can  be  defined. 

Real-Time  Processor  Architecture  Design 

A  real-time  implementation  architecture  for  the  algorithm 
needs  to  be  developed,  and  a  candidate  hardware  implementation 
should  be  identified.  Specifically,  the  following  issues  should 
be  addressed. 

•  Generation  of  an  architecture  design  that  best  meets 
the  features  of  the  detailed  algorithm  design  and  the 
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established  processor  requirements .  The  result  may  be 
an  architecture  with  features  different  from  those  in 
existing  processors,  and  which  is  likely  to  consist  of 
various  fundamental  architectures  (systolic;  vector; 
parallel  arrays;  etc.). 

•  Analysis  of  state-of-the-art  processors  to  determine 
which  contemporary  and  next-generation  VLSI  components 
best  match  the  optimized  architecture  design  and  the 
requirements . 

In  addressing  these  issues  the  emphasis  should  be  on  the  most 
computation-intensive  tasks  of  the  SSC  multichannel  algorithm 
reported  here . 

Processor  Development  System  Design 

A  processor  development  system  should  be  designed  and 
installed  to  serve  as  a  testbed  for  the  development  of  detection 
and  identification  methodologies  and  algorithms.  The  system 
should  be  applicable  for  on-line  laboratory  experimentation  and 
for  off-line  processing  of  data  collected  using  operational  radar 
systems.  Availability  of  such  a  development  system  will  speed  up 
significantly  the  algorithm  development  work  at  both  SSC  (during 
Phase  II)  and  RL  (after  delivery  upon  program  conclusion)  because 
the  generation  of  detection  results  require  Monte  Carlo  analyses 
and  simulations. 
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APPENDIX  A.  STATE  SPACE  REPRESENTATION  OP  TIME  SERIES 

MODELS 

Consider  a  discrete-time,  time-invariant,  complex-valued, 
zero-mean,  random  process  {&(n)}  defined  as  the  output  of  the 
following  state  space  system  model 

(A-ia)  y(n+1)»F¥(n)  +  Gu(n) 

(A-lbj  x(n)  *  HHy(n)  +  DHj^(n) 

Vector  recursive  processes  such  as  moving-average  (MA) ,  auto¬ 
regressive  (AR)  ,  and  auto-regressive  moving-average  (ARMA) 
processes  can  be  modeled  with  state  variable  models  (SVMs)  of  the 
form  (A-l) .  The  discussion  herein  is  limited  to  the  particular 
case  where  the  matrix  coefficients  of  the  recursion  are  square 
matrices,  and  the  number  of  output  coefficients  is  equal  to  the 
number  of  input  coefficients.  The  generation  of  a  minimal-order 
SVM  for  a  vector  recursive  process  involves  the  properties  of 
polynomial  matrix  pairs  and  canonical  forms  for  multiple  input, 
multiple  output  SVMs. 

In  contrast,  minimal-order  SVMs  for  scalar  recursive 
processes  (MA,  AR,  ARMA)  can  be  generated  in  a  straightforward 
manner  given  the  recursion  coefficients.  The  SVM  generic  form 
appropriate  for  modeling  scalar  recursive  processes  is 

(A-2a)  ^(n+1 )  -  Fy(n)  +  flu(n) 

(A-2b)  x(n)  *  hHx(n)  +  d*w(n) 

This  SVM  is  a  single-input,  single-output  system. 
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K .  1  Scalar  MX  Proenaa  Modal 


A  scalar  MA  process  of  order  M  is  defined  as 

x(n)  -  X  b*ku(n-k) 

k-0 

x(n)  ■  bju(n)  +  b*u(n-1)  +  bsu(n-2)  +  . . .  +  b^utn-M) 


where  {u(n)}  is  a  zero-mean  white  noise  sequence.  This  recursion 
can  be  modeled  with  a  state-space  system  of  the  form  (A-2)  with 
input  sequence  (u(n)},  and  state  vector  with  elements  that  are 
determined  by  the  input  sequence, 


'  yi(n)  " 

u(n-1) 

*(n)  - 

u(n-M+1) 

yu(n)  J 

u(n*M) 

The  output  noise  sequence  is  also  equal  to  the  input  noise 
sequence, 

w(n)  ■  u(n)  V  n 

which  means  that  the  input  and  output  noise  sequences  in  the  state 
space  model  are  completely  correlated.  Model  parameters  (F,  fl,  ft, 
d)  are  defined  as 


F 
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1 


a-ii- 


Qm-1 


HH-bH 

d*  •  b*0 

with  ^  denoting  a  IxM  vector  defined  by  M  of  the  MA  recursion 
coefficients, 

bH - [ b*  bj  ...  bj 

The  special  form  of  matrix  F  is  one  of  the  possible  four 
variations  of  the  so-called  companion  matrix  form.  Also,  the 
system  parameters,  the  quadruple  (F,  fl,  Jl,  d)  ,  is  a  variation  of 
the  so-called  controllable  canonical  form.  These  forms  have  the 
minimal  number  of  non-zero  elements  (whereby  the  name  "canonical") 
of  all  possible  SVMs  that  model  the  scalar  MA  process. 

Note  that  the  definition  of  the  state  vector  ^(n)  in  terms  of 
the  sequence  (u(n)}  inherently  defines  the  initial  condition  vector, 
*(0).  Once  the  initial  condition  vector  is  defined,  the  state 
propagation,  Equation  (A-2a) ,  provides  for  continued  generation  of 
the  output  process. 

Verification  of  the  above-defined  model  proceeds  as  follows. 
The  form  of  matrix  F  provides  for  continued  "scrolling"  of  the 
input  noise  sequence  as  elements  of  ^(n),  for  all  n.  Validation  of 
the  model  follows  from  (A-2b)  and  the  definition  of  ]l,  w(n),  and 
¥(n).  That  is, 
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x(n)  ■  JiHjt(n)  +  d*w(n)  =  b*V(n)  +  bju(n) 


Expanding  the  term  j2Hy(n),  and  substitution  of  the  definition  of  y(n) 
in  terms  of  the  sequence  {u(n)}  results  in 

x(n) «  bju(n)  +  biu(n-l)  +  b2u(n-2)  +  . . .  +  bMii(n-M) 

which  is  the  MA  process  definition.  Model  validation  can  be 
carried  out  also  using  the  transfer  function  concept,  as 
summarized  next . 

Consider  first  the  derivation  of  the  transfer  function  from 
the  MA  process  definition.  Since  the  MA  process  is  a  discrete¬ 
time  process,  the  appropriate  tool  for  the  determination  of  the 
transfer  function  is  the  Z-transform.  Application  of  the  Z- 
transform  to  the  definition  of  the  MA  model  results  in  the 
expression 

m  -  £  b^uu) 

k-0 


where  Z  denotes  the  transform  variable,  and  X(z)  and  U(z)  are  the  z- 
transforms  of  the  sequences  {x(n)}  and  {u(n)},  respectively.  The 
transfer  function  for  this  linear  system  is  then  defined  as 


T(z) 


X(z) 

U(z) 


M 

•XbU- 

k-0 


k 


This  corresponds  to  the  transfer  function  of  an  all-zero  system, 
as  is  well  known. 


The  transfer  function  for  a  single-input,  single-output  state 
variable  model  (A-2)  is  of  the  form 
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TM-tfkl-FTVd* 


The  particular  characteristics  of  matrix  F  and  vector  lead  to  a 
very  simple  expression  for  the  product  [zl  -  F]'1^;  namely, 


[zl  - 


m 


where  y(z)  is  the  system  characteristic  polynomial  (the  determinant 
of  matrix  [zl- F]) , 

7(2)  -  2M 

and  fi(z)  is  vector  with  elements  of  the  form  0j(z)  =  Z  ,*  that  is, 
fiT(z)  *  [  Z^1  •  •  •  Z2  Z  1  ] 

u  ^ 

Substitution  of  these  expressions  and  of  &  and  d  in  the  equation 
for  the  transfer  function  leads  to  the  following  result 

1(2),  b^M  .  z*M[bJ,zM  +  bHfl(z)] 

T(z)  -  b*o  +  b,  z’ 1  +  b‘2z*2  +  . . .  +  b'MzM  -  X  bkz'k 

k-0 


This  result  is  identical  to  the  transfer  function  expression 
derived  from  the  definition  of  the  MA  process. 

A.  2  Scalar  Alt  Pro e«n  Modal 

A  scalar  AR  process  of  order  M  is  defined  as 


9  1 


M 

x(n)  ■  *  X  a*kx(n-k)  +  u(n) 

k-1 

x(n)  ■  -  a^x(n-l)  -  a^x(n-2)  -  ...  -  aJ^n-M)  +  u(n) 

where  {u(n)J  is  a  zero-mean  white  noise  sequence.  This  recursion 
can  be  modeled  with  a  state-space  system  of  the  form  (A-2)  with 
input  sequence  {u(n)},  and  state  vector  with  elements  that  are 
determined  by  the  output  sequence, 


'  Yi(n)  ‘ 

x(n-1) 

X(n)  - 

yM>) 

m 

x(n-M+1) 

L  yM(n)  J 

x(n-M) 

The  output  noise  sequence  is  equal  to  the  input  noise  sequence, 

w(n)  ■  u(n)  V  n 


This  implies  complete  correlation  between  the  input 
noise  sequences  in  the  SVM  (as  in  the  case  of  the 
Model  parameters  (F,  fl,  ft,  d)  are  defined  as 


and  output 
MA  model) . 


92 


d*«  1 

with  &  denoting  a  vector  with  elements  equal  to  the  AR  recursion 
coefficients, 

aH-[a#i  ^  ...  a^] 

The  system  parameters  quadruple  (F,  g,  ft,  d)  ,  is  in  controllable 
canonical  form,  as  in  the  MA  model  case. 
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Note  that  the  definition  of  the  state  vector  Jt(n)  in  terms  of 
the  sequence  {x(n)}  inherently  defines  the  initial  condition  vector, 
m-  Once  the  initial  condition  vector  is  defined,  the  state 
propagation,  Equation  (A-2a) ,  provides  for  continued  generation  of 
the  output  process. 

Verification  of  the  above-defined  model  proceeds  as  follows. 
From  (A-2a)  and  the  definition  of  F,  y(n),  Jt(n+1),  and  fl,  it  follows 
that 


yM(n+1)  -  -  a*,yt(n)  -  a2y2(n)  -  ...  -  a*MyM(n)  +  u(n) 
yM(n+1)*-aH^(n)  +  u(n) 

Also,  it  follows  from  (A-2b)  and  the  definition  of  Jl,  w(n),  and  y(n) 
that 


x(n)  m  ftHy(n)  +  w(n)  » -  aHyi(n)  +  u(n) 

which  indicates  that  x(n)  ■  yM(n+1) .  Then,  expanding  the  term  •  aHX(n) 
and  substitution  of  the  definition  of  ^(n)  in  terms  of  the  sequence 
(x(n)>  results  in 

x(n)  ■  -  a^xfn-l)  *  a2x(n-2) «...  -  a]^x(n*M)  +  u(n) 

which  is  the  AR  process  definition. 

The  transfer  function  approach  can  be  used  also  to  validate 
this  SVM  for  scalar  AR  processes.  Application  of  the  Z-transform 
to  the  definition  of  the  AR  model  results  in  the  expression 
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M 

S 


a'krkX(z) .  U(z) 


where  a0  ■  1  is  introduced  for  notational  simplicity,  and  X(z)  and 
U(z)  are  the  Z-transforms  of  the  sequences  {x(n)}  and  {u(n)}, 
respectively.  The  transfer  function  for  this  linear  system  is 
then  defined  as 


T(z) 


m 

U(z) 


This  corresponds  to  the  transfer  function  of  an  all-pole  system, 
as  is  well  known. 


Consider  now  the  transfer  function  for  the  state  variable 
model  (A-2)  .  In  the  present  AR  process  case,  the  system 
characteristic  polynomial  is 

tfz)  -  zM  +  a'z^1  + . . .  +  a*M.,z  +  a’M 

and  the  particular  characteristics  of  matrix  F  and  vector  fl,  lead 
the  same  simple  expression  for  the  product  [zl  -  F]*1^;  namely, 

111 '  Fr’a  *  w  m 

where  &(z)  is  as  defined  previously.  Notice  that  the 
characteristic  polynomial  can  be  expressed  as 

y(z)  -  zM  +  aHfi(z) 

Substitution  of  these  expressions  and  of  and  d  in  the  equation 
for  the  transfer  function  leads  to  the  following  result 
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T(z) 


hHfl(z)  +  d*y(z)  _^fi(z)jM((z[  2m 
7(z)  "  y(z)  “  y(z) 


T(z) 


1 

2*M7(2) 


This  is  identical  to  the  transfer  function  expression  derived  from 
the  definition  of  the  AR  process . 

a. 3  aeiAix  ABMA _ Laataa.  Modal 

A  scalar  ARMA  process  of  order  M  is  defined  as 

M  M 

x(n)  -  *  X  a*kx(n-k)  +  X  t>ku(n*k) 

k-i  k-o 

x(n)  -  *  a^xfn-l)  -  ...  -  aJ/l.1x(n-M+1)  -  a^x(n-M)  +  boii(n)  +  b^u(n-1 )  + 

+  b2U(n-2)  +  . . .  +  b^uln-M) 

where  {u(n)}  is  a  zero-mean  white  noise  sequence.  This  recursion 
can  be  modeled  with  a  state-space  system  of  the  form  (2)  with 
input  sequence  {u(n)},  and  output  noise  sequence  equal  to  the  input 
sequence, 


w(n)  *  u(n)  V  n 

This  implies  complete  correlation  between  the  input  and  output 
noise  sequences  in  the  SVM  (as  in  the  case  of  the  MA  and  the  AR 
models)  .  Model  parameters  (F,  g,  ft,  d)  are  defined  as 
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r  i  i 


HH  =  bH- b*oaH 

d*  -  b*o 

u 

Here,  as  in  the  AR  case,  vector  a  has  elements  equal  to  the  AR 
recursion  coefficients, 

SH  *  [  a1  *2  1  •  •  ^  ] 

m 

and  vector  b.  has  elements  defined  by  M  of  the  MA  recursion 
coefficients, 
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bH  *  [  bi  bj  ...  bM] 


The  system  parameters  quadruple  (F,  Q.,  Jl,  d) ,  is  in  controllable 
canonical  form,  as  in  the  MA  and  AR  model  cases. 

State  vector  initial  conditions,  y,(0),  for  this  case  are 
related  to  the  input  and  output  sequences  in  a  more  complex 
manner,  and  have  to  be  selected  appropriately.  Once  the  initial 
condition  vector  is  defined,  the  state  propagation,  Equation  (A- 
2a),  provides  for  continued  generation  of  the  output  process. 


The  simplest  approach  to  validate  this  model  is  via  the 
transfer  function  approach.  Application  of  the  Z-transform  to  the 
definition  of  the  ARMA  model  results  in  the  expression 

£  a'krkx(z)  -  £  bkz-kU(z) 

k-0  ksO 


where,  as  before,  X(z)  and  U(z)  are  the  Z-transforms  of  the 
sequences  {x(n)}  and  {u(n)},  respectively,  and  Sq  ■  1  is  introduced  for 

notational  simplicity.  The  transfer  function  for  this  linear 
system  is  then  defined  as 


M  M 

y  bU*k  y  b!zM'k 

T(z)  ■  X<z>  .  k-°  -k-° 

U(Z)  M  M 

Z  <2‘k  £  a;zM’k 

k-0  k-0 


where  the  two  polynomial  ratio  expressions  (corresponding  to 
inverse  powers  of  Z  or  direct  powers  of  Z)  are  equivalent,  as 
indicated.  This  is  a  transfer  function  with  both  poles  and  zeros, 
as  expected  for  an  ARMA  process. 
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Consider  now  the  transfer  function  for  the  state  variable 
model  (A-2) .  For  an  ARMA  process  the  system  characteristic 
polynomial  is 

Y(z)  -  zM  +  a;  z^1  + . . .  +  z  +  a*M 

which  is  equal  to  that  for  an  AR  process  SVM  model.  As  in  the 
other  two  cases, 

(zl  •  .  -L-  6(2) 

given  the  particular  features  of  matrix  F  and  vector  fl.  (fi(z)  is  as 
defined  previously) .  Notice  also  that,  as  in  the  AR  process  case, 
the  characteristic  polynomial  can  be  expressed  as 

tfz)-  zM  +  aHfi(z) 

it  * 

Substitution  of  these  expressions  and  of  tl  and  d  in  the  equation 
for  the  transfer  function  leads  to  the  following  result 

T/  v  llHfl(z)  d*Y(z)  (bH  -  b*oaH)  ft(z)  +  b*o7(z)  bpzN  +  bHQ(z) 

V  Y(z)  "  .  Y(z)  "  Y (z) 

It  is  easy  to  verify  that  this  result  is  identical  to  the  transfer 
function  expression  derived  from  the  definition  of  the  ARMA 
process.  That  is, 


M 

y 

tm,  b'nz”  ■»  hHfi(z)  t,0  k 

Y(z)  m 

I  <z”* 

k-0 
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where  SLq ■  1 ,  as  before. 


X.  4  Modala  for  Vector  Reeuraive  Proctaasa 

Vector  recursive  processes  of  the  MA,  AR,  and  ARMA  type  can 
be  represented  with  SVMs  of  the  type  given  herein.  For  vector 
recursive  processes  the  appropriate  notation  is: 

ma  X(n)  ■  Bju(n)  +  B”u(n-1 )  +  B^n-a)  +  . . .  +  B[Ju(n-M) 

ar  X(n)  -  -  A”x(n-1 )  -  Aji(n-2)  -  ...  *  Afc(n-M)  +  u(n) 

arma  *(n)  -  -  A*fe(n-1)  -  . . .  •  A^^n-M+I)  *  A^n-M)  +  Bfti(n)  +  B?u(n-1 )  + 

+  Bgu(n-2)  +  . . .  +  BjJu(n-M) 

where  each  of  the  coefficient  matrices  is  dimensioned  JxJ .  Also 
analogous  to  the  scalar  case,  the  corresponding  transfer  function 
matrices  can  be  defined  using  the  z-transform;  which  leads  to 

Tma(z)  -  B(z) 

Tar(z)«A-’(z) 

WW-a-’wbw 

where  A(z)  and  B(z)  are  the  following  matrix  polynomials  in  Z, 

A(z)«  £  A?z-k 

K-0 

B(2)-I  Bfz-k 

k-0 
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with  Aq  the  JxJ  identity  matrix.  The  matrix  pair  (A(z),  B(z)} 
(including  the  cases  with  either  A(z)  ■  I  or  B(z)-I)  corresponding  to 
a  linear  discrete-time  system  is  referred  to  as  a  matrix 
polynomial  representation  or  a  matrix  fraction  description  (MFD) 
for  the  system. 

Departing  from  the  time-domain  definition  for  the  vector 
recursive  processes,  the  SVM  for  each  of  the  three  processes  is  of 
the  same  form  as  the  corresponding  scalar  case  SVM,  with  the 
following  changes:  a  JxJ  coefficient  matrix  in  place  of  the 
corresponding  coefficient  scalar,  a  JxJ  identity  matrix  (lj)  in 
place  of  each  unit  scalar,  and  a  JxJ  null  matrix  (Oj)  in  place  of 

each  zero-valued  scalar.  Specifically,  the  SVM  for  the  ARMA 
vector  process  is : 
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L  o,  J 

HH  -  [  B?  ■  A?B”  B»  ■  AjBq 
Dh  ■  Bq 


The  SVM  for  the  other  vector  processes  (MA;  AR)  is  obtained  by 
substituting  the  correct  values  for  the  vector  process 
coefficients  in  the  above  system  parameters  (that  is,  A| » Oj  for  an 
MA  process;  and  Bq  »  lj  and  B|  ■  Oj,  i  £  1  for  an  AR  process).  In  all 
cases,  the  transfer  function  matrix  is  obtained  from  the  SVM 
representation  as 

T(z)  -  Hh[zI  -  F^G  +  Dh 


A  transfer  function  calculated  according  to  this  relation  is 
equivalent  to  the  transfer  function  calculated  from  the 
appropriate  polynomial  matrices. 

The  order  (dimension  of  the  state  vector)  of  the  resulting 
SVM  for  each  of  the  three  vector  processes  is  N  *  MJ,  since  for 
each  process  the  system  matrix  F  consists  of  M  block  rows  and  M 
block  columns,  where  each  block  in  each  row  and  column  is  a  JxJ 
matrix.  SVM  order  is  important  for  practical  and  computational 
considerations.  An  SVM  representation  is  of  minimal  order  if  no 
other  SVM  representation  of  lower  order  leads  to  the  same  transfer 
function  matrix.  In  terms  of  the  system  parameters  (F,  G,  H,  D) , 
the  order  of  the  SVM  representation  is  determined  by  the  rank  of 
the  controllability  matrix  or  the  rank  of  the  observability 
matrix,  whichever  is  smaller.  Given  the  form  of  the  matrix  pair 
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(F.G)  for  all  three  cases,  it  is  easy  to  verify  that  the 
controllability  matrix  has  full  rank  for  all  three  cases. 
However,  the  observability  matrix  has  a  simple  form  only  for  the 
MA  SVM.  The  special  form  of  the  observability  matrix  for  the  MA 
case  considered  herein  (with  a  square  matrix)  indicates  by 

inspection  that  the  rank  of  the  observability  matrix  is  equal  to 
MJ  if  and  only  if  matrix  has  full  rank.  Such  a  simple  result 

is  not  available  for  the  AR  and  the  ARMA  SVMs.  Determination  of 
the  conditions  on  the  coefficients  of  the  polynomial  matrices  A(z) 
and  B(z)  for  AR  and  ARMA  vector  processes  that  lead  to  an  SVM 
representation  of  minimal  order  is  a  difficult  problem.  This  is 
due  to  the  fact  that  both  AR  and  the  ARMA  vector  processes  lead  to 
a  transfer  function  matrix  with  elements  which  are,  in  general,  a 
ratio  of  polynomials . in  Z. 

Model  order  and  related  issues  for  matrix  polynomial 
representations  have  been  discussed  by  several  researchers.  The 
results  summarized  next  are  available  in  the  text  by  Rosenbrock 
(1970).  Consider  the  matrix  polynomial  representation  of  a 
system,  and  assume  that  the  determinant  of  A(z)  is  different  from 
zero  to  eliminate  pathological  cases.  For  an  AR  vector  process, 
the  order  of  the  system  is  given  by  the  degree  of  the  determinant 
of  A(z) .  Thus,  the  SVM  representation  presented  herein  for  vector 
AR  processes  is  of  minimal  order  if  the  determinant  of  A(z)  (with 
A0  ■  lj)  has  degree  equal  to  MJ. 

Several  definitions  need  to  be  introduced  prior  to  stating 
the  relevant  results  regarding  minimal  order  for  ARMA  vector 
processes.  A  square  polynomial  matrix  is  said  to  be  regular  when 
the  matrix  coefficient  of  the  highest  power  of  Z  is  non-singular. 
The  determinant  of  a  regular  polynomial  matrix  has  maximum 
possible  degree.  A  square  polynomial  matrix  is  said  to  be 
unlmodular  if  its  determinant  is  a  non-zero  constant.  Unimodular 
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polynomial  matrices  have  an  inverse  which  is  also  a  polynomial 
matrix.  As  an  example,  the  polynomial  matrix 


Q(z)  *  Q0  +  C^z'1 


1  +  z'1  3  +  z'1 

2  +  z*1  4  +  z'1  . 


is  unimodular  because  the  determinant  of  Q(z)  is  equal  to  -2. 
Notice  that  the  inverse  of  Q(z)  is  also  a  polynomial  matrix, 


0-’(z)--L 


4  +  z'1 
_  -  (2  +  z'1 ) 


-  (3  +  z'1 ) 

1  +Z1 


as  expected.  Notice  also  that  Q(z)  is  not  a  regular  matrix  since 
Ql  is  singular. 


Two  polynomial  matrices  A(z)  and  B(z)  are  said  to  have  a  common 
left  divisor  $(z)  if 

A(z)  =  S(z)PA(z) 

B(z)  =  S(z)PB(z) 

where  S(z),  PA(z),  and  PB(z)  are  polynomial  matrices.  Finally,  if 
all  the  common  (left)  divisors  of  two  polynomial  matrices  A(z)  and 
B(z)  are  unimodular,  then  the  two  matrices  are  said  co  be 
relatively  (left)  prime.  That  is,  if  A(z)  and  B(z)  are  relatively 
(left)  prime,  then  the  determinant  of  the  polynomial  matrix  S(z)  in 

the  above  factorizations  is  a  constant.  This  implies  that  the 
degree  of  the  determinant  of  PA(z)  is  equal  to  the  degree  of  the 
determinant  of  A(z),  and  the  degree  of  the  determinant  of  PB(z)  is 
equal  to  the  degree  of  the  determinant  of  B(z) .  Furthermore,  the 
determinant  of  PaU)  has  no  polynomial  factors  in  common  with  the 
determinant  of  PB(z).  A  matrix  polynomial  pair  (A(z),  B(z))  with  A(z) 
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and  B(z)  relatively  (left)  prime  is  an  irreducible  matrix 
polynomial  representation  for  the  system. 

The  relevant  results  for  ARMA  vector  processes  can  be  stated 
now.  As  in  the  AR  case,  for  an  ARMA  vector  process  the 
determinant  of  A(z)  (with  A0  a  lj)  must  have  degree  equal  to  MJ  for 

the  SVM  representation  presented  herein  to  be  of  minimal  order. 
However,  two  additional  conditions  must  be  satisfied.  Namely, 
matrix  must  have  full  rank,  and  the  polynomial  matrices  A(z)  and 
B(z)  must  be  relatively  (left)  prime.  Full  rank  for  matrix  B^ 
implies  that  B(z)  is  a  regular  polynomial  matrix.  If  A(z)  and  B(z) 
are  not  relatively  prime,  then  the  order  of  the  system  is  reduced 
by  the  degree  of  the  determinant  of  the  greatest  common  (left) 
divisor  of  A(z)  and  B(z) .  This  is  related  to  the  so-called 
pole/zero  cancelations. 
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APPENDIX  B.  QUOTIENT  SINGULAR  VALUE  DECOMPOSITION 


The  quotient  singular  value  decomposition  (QSVD)  is  a 
generalization  of  the  SVD  for  a  matrix  to  the  case  of  two  general 
matrices  with  the  same  number  of  columns.  As  such,  it  is  also 
referred  to  as  the  generalized  SVD.  This  concept  was  developed  by 
Van  Loan  (1976),  who  called  it  the  BSVD .  Paige  and  Saunders 
(1981)  modified  the  concept  and  extended  its  applicability  to 
general  matrices.  Their  concept  is  summarized  herein  in  the 
context  of  the  multichannel  detection  application.  From  a 
computational  viewpoint,  the  approach  suggested  by  Van  Overschee 
and  De  Moor  (1993)  is  adopted.  Two  distinct  cases  are  considered 
herein,  corresponding  to  the  two  conditions  that  arise  in  the 
implementation  of  the  identification  algorithm  (Section  3.1). 


B.  1  QSVD _ for _ the  Matrices  of  Equations  (3-22)  and  (3-231 

Consider  the  J(L-1)xJ(L+1)  matrix  RD,  and  the  J(L-1)xJ(L-1)  matrix 
Re  defined  in  Equations  (3-22)  and  (3-23),  respectively.  It  is 

desired  to  determine  the  QSVD  of  the  matrix  pair  consisting  of  the 
conjugate  transpose  of  these  two  matrices.  The  procedure  is 
described  below. 


The  first  step  is  to  define  a  2JLxJ(L-1)  matrix  A  as  the 
following  concatenation  of  the  conjugate  transposes  of  matrices  R0 
and  Re: 


(B-l) 


A  = 


Now  carry  out  an  SVD  on  matrix  A  to  get  (recall  that  A  has  more 
rows  than  columns) , 
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<B-2a) 


A  -  U 


SA 

101  . 


UlA 

U2A 


U3A 

U4A 


u 

u 


- 

'S1A 

[0]  ' 

5A 

[0] 

[0] 

6A  - 

.  [0] 

[0]  . 

(B-2b) 


S,A  [01 
.  [01  [0]  . 


In  this  decomposition  the  unitary  matrix  UA  is  J(L+1 )xJ(L+1 ),  matrix 
SA  is  J(L-1)xJ(L*1),  and  the  unitary  matrix  VA  is  J(L-1  )xJ(L*1 ) .  Matrix 
SA  is  diagonal,  with  real-valued  non-negative  elements  along  the 
diagonal  arranged  in  decreasing  order  of  magnitude  (the  largest- 

valued  element  occupies  the  (1,1)  position)  .  The  diagonal 
elements  of  matrix  SA  are  the  singular  values  of  matrix  A.  The 

rank  of  matrix  A,  denoted  herein  as  KA  =  rank (A),  is  equal  to  the 

number  of  non-zero  singular  values.  These  non-zero  singular 
values  are  the  diagonal  elements  of  the  KAxKA  matrix  S1A.  If 

matrix  A  is  full-rank,  then  S1A  becomes  SA.  In  most  cases 
involving  random  processes  the  singular  values  of  A  will  be  non¬ 
zero,  although  there  may  be  a  large  dynamic  range  between  the 
largest  and  the  smallest  singular  values.  Alternatively,  the 
singular  values  may  appear  in  groups,  with  a  significant  variation 

in  dynamic  range  between  the  groups  of  singular  values.  Such 
situations  are  indicative  of  an  effective  rank  KA  <  J(L-1). 

The  row  partition  of  matrix  UA  into  submatrices  is  the  same 

as  the  row  partition  in  the  concatenation  (B-l),  and  the  column 
partitioning  of  UA  corresponds  to  the  row  partitioning  of  matrix 

SA.  That  is,  matrix  U1A  is  J(L+1)xKA  and  matrix  U2A  is  J(L-1)xKA. 

With  this  partitioning,  Equation  (B-2a)  is  equivalent  to  the 
following  expression, 
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(B-3)  A  - 

■u1A  NT 

'  s1A  10]  ' 

VA  * 

■U,A  [0]  • 

- 1 

CO 

> 

o 

f 

.Uza  [01. 

.  [0]  [0]  . 

,U2A  [0], 

L  [01  IM  J 

where  M*J(L-1)-KA.  This  expression  is  obtained  by  first  cutting 
off  the  last  block  column  of  matrix  UA  and  the  corresponding  block 
row  of  matrix  SA.  Next,  U3A  and  U4A  are  replaced  with  zeros,  which 
allows  placing  an  identity  matrix  of  dimensions  M»J(L-1)-KA  into 
the  lower  right-hand-corner  of  matrix  SA.  None  of  these 

modifications  alters  the  numerical  value  of  the  expression. 

The  next  step  in  the  computation  of  the  QSVD  is  to  carry  out 
an  SVD  on  each  of  the  matrices  U1A  and  U2A.  The  resulting  SVDs  can 

be  expressed  as 

(B-4)  U1A  »  UL..,  Sl..,  Vy1A 

(B-5)  »  VL>1  Tw  Vu2a 

In  the  first  decomposition,  UL.i  is  a  J(L+1  )xJ(L+1 )  unitary  matrix, 
SL..,  is  J(L+1)xKa,  and  Vy1A  is  a  KAxKA  unitary  matrix.  Matrix  SL.i  is 

zero  except  for  real-valued,  non-negative  elements  along  the  main 
diagonal.  The  non-negative  elements  of  S|_„i  are  arranged  in 

decreasing  order  of  magnitude,  with  the  largest-valued  element 

occupying  the  (1,1)  position  and  having  value  less  than  or  equal 
to  unity.  In  the  second  decomposition,  is  a  J(L*1  )xJ(L.-1 ) 

unitary  matrix,  TL..,  is  J(L-1)xKA,  and  VU2A  is  a  KAxKA  unitary 
matrix.  Matrix  is  zero  except  for  real-valued,  non-negative 

elements  along  the  main  diagonal.  The  non-negative  elements  of  T|_m1 

are  arranged  in  increasing  order  of  magnitude,  with  the  smallest- 

valued  element  occupying  the  (1,1)  position  and  having  value 
greater  than  or  equal  to  zero.  The  largest-valued  element  of 

is  the  (Ka, Ka) th  element,  and  its  value  is  less  than  or  equal  to 
unity.  Notice  that  the  arrangement  of  the  elements  of  T«  1  along 
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the  main  diagonal  is  the  reverse  of  the  conventional  SVD.  This 

deviation,  however,  is  essential  to  the  Paige  and  Saunders  (1981) 
definition  of  the  QSVD.  Let  {Si  |  i  ■  1,  2,  ...  ,  KA}  denote  the  main 

diagonal  elements  of  SL.i,  and  let  {tj  |  i  »  1 , 2,  ...  ,  KA)}  denote  the  main 
diagonal  elements  of  TL.., .  The  above-stated  conditions  on  these 

elements  are  summarized  as: 

(B-6)  1  ^S1  2S2£. . .  £SKa£0 

(B-7)  0  S  t1  S  t2  S  .  .  .  S  tKA  Si 

Paige  and  Saunders  (1981)  ha/e  shown  also  that  these  elements 
satisfy  the  following  constraint, 

(B-8)  sf  +  tf-1  i  =  1 , 2,  . . .  ,  Ka 

This  constraint  is  valid  only  if  the  singular  values  satisfy  (B-6) 
and  (B-7).  The  pairs  of  values  (Sj,  tj)  are  called  the  singular  value 
pairs  of  matrices  Rg  and  RE. 

Based  on  Equations  (B-6)  and  (B-7)  and  on  the  orthogonality 
property  of  matrix  UA  it  is  possible  to  show  that 

(B-9)  VU1A>VU2A  =  VUA 


Then,  substituting  this  result  into  Equations  (B-4)  and  (B-5),  and 
in  turn  substituting  these  into  Equation  (B-3)  leads  to 


(B-10) 

A  = 

[0]  ■ 

- 1 

< 

c 

> 

[0]  ' 

'S1A  [0]  ' 

[0]  . 

L  [0] 

- 

.  [0]  |M  . 

Now  define  a  J(L-1  )xJ(L-1 )  matrix  Yl.i  as 
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(B-ll) 


3 

> 

i _ 

[0]  ' 

‘S1A 

[0]  ' 

[0] 

.  10] 

>M  . 

Substitute  for  Yl.i  in  Equation  (B-10)  and  re-arrange  the 
submatrices  of  the  first  matrix  to  obtain 


(B-12)  A 


uL.t  [  sL.t 

VL.,  [  Tl, 


The  desired  QSVD  for  the  matrix  pair  Rq  and  Rg  follows  directly  by 
a  comparison  of  Equations  (B-l)  and  (B-12), 

(B-i3)  Ro-ul.,[sl.,  oJIU,)M]y^, 

(B-14 )  Rg  »  VL., [ Tl.,  Oj(L.,1m ]  Y», 

If  matrix  A  has  full  rank  these  expressions  simplify  to 
(B-15)  Rg  -  Uu.,SL.,Yg, 

(B-l  6)  Rg  -  VL.,TL.,Yg., 


As  mentioned  earlier,  in  most  cases  involving  random  processes 
matrix  A  will  be  full  rank.  Even  if  such  is  not  the  case,  it 
appears  to  be  better  to  over-estimate  the  rank  of  A  rather  than  to 
under-estimate  it.  In  fact,  if  the  rank  of  matrix  A  is  under¬ 
estimated,  then  the  true  structure  of  matrices  and  is 

distorted.  With  over-estimation  of  the  rank  of  A  it  is  still 
possible  to  determine  the  true  model  order  accurately,  but  under¬ 
estimation  of  the  rank  of  A  effectively  places  an  upper  bound  on 
the  attainable  model  order  and  this  bound  could  be  less  than  the 
true  model  order  (see  Section  3.2) . 


no 


Consider  matrix  for  the  general  case,  where  rank(A)  ■  KA  S 

J(L*1).  The  structure  of  this  matrix  is  determined  by  the  true 
order  of  a  state  space  representation  of  the  process  being 
modeled.  Specifically,  for  the  case  where  the  model  order  is  N  < 
KA'  matrix  SL.i  is  of  the  form 


rod) 

bL-1 

[0]  “ 

~  ed) 
bL-1 

[of 

(B-17) 

SL.i  - 

[0] 

e(2> 

bL-1 

. 

[0] 

[0] 

.  [0] 

[0]  . 

.  [0] 

[0]. 

Here  S[^  is  an  NxN  diagonal  matrix,  and  S[2^  is  a  (KA-N)x(KA-N) 

matrix  with  possible  non-zero  elements  only  along  the  main 

(2) 

diagonal.  As  inferred  by  Equation  (B-17)  ,  is  a  null  matrix 

when  the  model  order  N  <  KA,  In  practical  situations  where 

(2) 

randomness  is  present,  the  diagonal  elements  of  matrix  are  not 

equal  to  zero,  but  they  are  significantly  smaller  than  the 
diagonal  elements  of  S^1 1 . 

B .  2  QSVD _ for  the  Matrices  of  Equations  (3-201  and  (3-211 

Consider  the  JLxJL  matrix  RB,  and  the  JLxJL  matrix  Rc  defined 

in  Equations  (3-20)  and  (3-21),  respectively.  It  is  desired  to 
determine  the  QSVD  of  the  matrix  pair  consisting  of  the  conjugate 
transpose  of  these  two  matrices .  Since  the  approach  is  analogous 
to  the  preceding  section,  only  the  key  steps  and  definitions  are 
given  below. 

As  before,  define  a  2JLxJL  matrix  B  as  the  following 
concatenation  of  the  conjugate  transposes  of  matrices  R0  and  Rc: 


(B-18) 


B 


'  «e' 
-  Rc. 


Now  carry  out  an  SVD  on  matrix  B  to  get  (recall  that  B  has  more 
rows  than  columns) , 


(B-19a) 


B  -  Uf 


»B 


[0]  J 


/H 


UlB 

U2B 


U3B 

u4B 


U5B 

U6B- 


Sib 

[0] 

.  [0] 


[0] 

[0] 

[0] 


j 


(B-19b) 


Sib  [0] 

.  [0]  [0]  . 


In  this  decomposition  the  unitary  matrix  Ug  is  2JLx2JL,  matrix  Sg 
is  2JLxJL,  and  the  unitary  matrix  Vg  is  JLxJL.  Matrix  Sg  is 

diagonal,  with  real-valued  non-negative  elements  along  the 
diagonal  arranged  in  decreasing  order  of  magnitude  (the  largest- 

valued  element  occupies  the  (1,1)  position).  The  diagonal 
elements  of  matrix  Sg  are  the  singular  values  of  matrix  B,  and  the 

rank  of  matrix  B,  denoted  herein  as  Kg  =  rank(B),  is  equal  to  the 

number  of  non-zero  singular  values.  These  non-zero  singular 
values  are  the  diagonal  elements  of  the  KgxKg  matrix  S1B.  If 

matrix  B  is  full-rank,  then  S1B  becomes  SB.  Analogous  to  the 

prior  case,  for  random  processes  the  effective  rank  of  matrix  B  is 

Kb  £  JL. 


As  before,  Equation  (B-19a)  can  be  converted  to  the  following 
equivalent  form, 


(B-20)  B  « 

*  U1B  [0]  * 

'  $,B  [0]  ' 

|0]  1 

1 

'S,B  [0]  ■ 

-  U2B  [0]  . 

.  10]  [0]  j 

D 

1 

a 

to 

CD 

o 

.  [0]  lM  . 
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(B-20)  B  * 

‘  U1B  [0]  ' 

'  S1B  [0]  ' 

VB  = 

■  U,B  [0]  ‘ 

‘s10  [0]  " 

o 

CD 

CM 

z> 

.  [0]  [0]  . 

.  U2B  [01  . 

.  [0]  >M  - 

where  M*JL-K0.  The  SVD  of  the  JLxK0  matrix  U10  and  the  SVD  of  the 
JLxK0  matrix  U20  are  expressed  as 

(B-21)  U10  ■  UL  SLVJ10 

(B-22)  U20  «  ^L^”l^U2B 

In  the  first  decomposition,  UL  is  a  JLxJL  unitary  matrix,  SL  is 
JLxK0,  and  Vjj10  is  a  K0xK0  unitary  matrix.  Matrix  SL  is  zero 

except  for  real-valued,  non-negative  elements  along  the  main 
diagonal.  The  non-negative  elements  of  SL  are  arranged  in 

decreasing  order  of  magnitude,  with  the  largest-valued  element 

occupying  the  (1,1)  position  and  having  value  less  than  or  equal 
to  unity.  In  the  second  decomposition,  VL  is  a  JLxJL  unitary 

matrix,  Tj_  is  JLxK0,  and  V^b  a  ^0xK0  unitary  matrix.  Matrix  T|_ 

is  zero  except  for  real-valued,  non-negative  elements  along  the 
main  diagonal.  The  non-negative  elements  of  T(_  are  arranged  in 

increasing  order  of  magnitude,  with  the  smallest-valued  element 

occupying  the  (1,1)  position  and  having  value  greater  than  or 
equal  to  zero.  The  largest-valued  element  of  T|_  is  the  (K0,K0)th 

element,  and  its  value  is  less  than  or  equal  to  unity.  As  in  the 
prior  case,  the  arrangement  of  the  elements  of  along  the  main 

diagonal  is  essential  to  this  definition  of  the  QSVD .  The  non¬ 
zero  elements  of  and  Tj_  are  the  singular  value  pairs  of  matrices 

R0  and  Rc,  and  they  satisfy  conditions  identical  to  (B-6)  through 
(B-8)  with  Ka  replaced  by  K0. 


Given  the  conditions  satisfied  by  the  singular  value  pairs 
and  given  the  orthogonality  property  of  matrix  U0  it  is  possible 

to  show  that 
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(B-23) 


VU1B  *  VU2B  =  VUB 


Then,  substituting  this  result  into  Equations  (B-21)  and  (B-22), 
and  in  turn  substituting  these  into  Equation  (B-20)  leads  to 


(B-24) 


B  - 

•  ulsl 

[01  ■ 

■vUB 

[01  ' 

Sib 

[01  ' 

.  VJl 

[0]  . 

.  [01 

lM . 

.  [01 

'm  . 

Now  define  a  JLxJL  matrix  YL  as 


(B-25) 

yl  - 

< 

C 

CD 

[01  ' 

Sib 

[0]  ' 

[0] 

•m  . 

.  [01 

'm  . 

Substitute  for  YL  in  Equation  (B-24)  and  re-arrange  the  submatrices 
of  the  first  matrix  to  obtain 


(B-26) 


'  UJ  SL 

.  vl[tl 


vr 


The  desired  QSVD  for  the  matrix  pair  Rg  and  Rc  follows  directly  by 
a  comparison  of  Equations  (B-18)  and  (B-26), 

(B-27)  R^  ■  UL[SL  0jlm]y£ 

(B-28)  RS-VJT,.  Ojl.m ] Yf 


If  matrix  B  has  full  rank,  as  can  be  expected  in  most  cases  where 
the  processes  are  random,  these  expressions  simplify  to 

(B-29)  Rg  -  Ul$lYl 


1  14 


Tht  commtnts  msdt  tirli«k  regarding  tha  *'*Nk  »' i  mat i In  A  apply 
•qually  to  tha  rank  of  matrix  B,  In  particular,  It  la  Tartar  t* 
ovar-aatimata  tha  rank  of  B  rathar  than  to  undar-aatimata  it, 


Considar  matrix  for  the  casts  whara  rank(B)  ■  Kg  £  JL,  Tha 

structure  of  thia  matrix  ia  datarminad  by  tha  trua  order  of  a 
stata  space  rapraaantation  of  the  process  being  modeled,  in  tact, 
modal  order  can  be  datarminad  by  examining  tha  diagonal  elements 
of  SL  (sea  Section  3.2).  Modal  order  can  be  datarminad  also  from 

tha  diagonal  elements  of  matrix  SL,, .  However,  it  is  preferable  no 
use  matrix  S[_  for  model  order  determination  because  this  matrix  is 

generated  by  the  QSVD  of  two  matrices  with  JL  rows.  i,«ca  a  is 

more  robust  numerically  than  the  QSVD  for  S^»  which  is  a  Q3VD  for 

two  matrices  with  J(L*1)  rows.  For  the  cases  where  the  modei  order 
is  N  <  Kg,  matrix  SL  is  of  the  form 


's<1) 

[0]  ' 

tor 

(B-31) 

SL  - 

[0] 

e(2) 

m 

[0] 

[0] 

-  [0] 

[0]  . 

.  [0] 

[0]. 

Here  is  an  NxN  diagonal  matrix,  and  S[2>  is  a  (Kg-N)x(Kg-N) 

matrix  with  possible  non-zero  elements  only  along  the  main 

(2) 

diagonal.  As  inferred  by  Equation  (B-31),  is  a  null  matrix 

when  the  model  order  N  <  KQ .  When  random  processes  are  being 

(2) 

modeled  the  diagonal  elements  of  matrix  are  not  equal  to  zero, 

but  their  magnitude  is  smaller  than  the  diagonal  elements  of  Sj^. 

In  such  cases  the  relative  numerical  value  of  the  elements  along 
the  main  diagonal  of  determines  the  cut-off  point,  and 

consequently,  the  model  order  (Section  3.2)  . 
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